# Gravity data inversion of the pyrenees range using taguchi sensitivity analysis and ADMM bound constraints based on seismic data

Roland Martin, Vitaliy Ogarko, Jérémie Giraud, Bastien Plazolles, Paul Angrand, Sonia Rousse, Melina Macouin

► To cite this version:

Roland Martin, Vitaliy Ogarko, Jérémie Giraud, Bastien Plazolles, Paul Angrand, et al.. Gravity data inversion of the pyrenees range using taguchi sensitivity analysis and ADMM bound constraints based on seismic data. Geophysical Journal International, 2024, 10.1093/gji/ggae410. insu-04831460

HAL Id: insu-04831460

https://insu.hal.science/insu-04831460v1

Submitted on 11 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Gravity data inversion of the Pyrenees range using Taguchi sensitivity analysis and ADMM bound constraints based on seismic data

Roland Martin,<fnref n="1" /> Vitaliy Ogarko,<fnref n="2" />,<fnref n="3" /> Jérémie Giraud,<fnref n="2" />,<fnref n="4" /> Bastien Plazolles,<fnref n="1" /> Paul Angrand,<fnref n="1" /> Sonia Rousse<fnref n="1" /> and Melina Macouin<fnref n="1" />

Accepted 2024 November 7. Received 2024 November 5; in original form 2023 December 13

## SUMMARY

Gravity inversion methods are able to recover density distributions in the Earth but they require strong constraints based on various prior information. In the present study, our aim is to invert gravity anomalies using existing geological models and density information in orogenic areas like the extensively studied Pyrenees region. Geophysical exploration for resource assessment using gravity and seismic data already accessible in this area has been performed to identify potential fluid/gas reservoirs of economic significance. For a given pre-existing model defined by different geological units and summarizing our knowledge of the area, we aim at retrieving the density distribution within each unit. For this, we use an Alternate Direction Multiplier Method to perform gravity data inversion by constrained by Interval Bound Constraints (IBC-ADMM) defined as bounding distribution intervals of possible density values. To estimate these bounding intervals, we first use a prior density model geologically compatible and obtained using geological information, gravity modelling, seismological data and seismic models. In a second step, we apply a Taguchi statistical analysis on representative density variations inside each unit of the prior density model to estimate their impact on data residuals, and reduce the parameter space prior to inversion. We perform the gravity data inversion constrained by these Taguchi-derived density intervals and we estimate model uncertainties for these density intervals. We apply the technique to the entire Pyrenees range at 2 km resolution. The sensitivity analysis shows that the strongest variations on gravity data are obtained essentially in the Axial Zone, in the middle and lower crust, and in the Upper mantle. Inverted densities are compared to the prior model. They are increased in the central part of the Pyrenees Axial Zone while, in the deeper structures (from middle crust down to upper mantle), they are decreased in the Southern (Spanish) Iberian upper mantle and increased in the Northern (French) Eurasian Upper mantle. A possible shortening of the Iberian lower crust–upper mantle system can be assumed below the western part of the Axial Zone.

**Key words:** Gravity anomalies and Earth structure; Inverse theory; Numerical modelling; Statistical methods; Wavelet transform; Crustal structure.

## 1 INTRODUCTION

Geophysical data inversion has long been a means to image the subsurface and complement geological models. In the last decade, the integration of both geophysical information and geophysical inversion has been focusing a lot of interest. Taking geological information in geophysical inversions allows to regularize the inverse

problem and obtain models physically and geologically compatible and consistent. To reduce model uncertainties, different approaches have been developed to incorporate relevant geological information, correlations between parameters or combinations of different physics. Furthermore, petrophysical information can also be introduced through clustering techniques. To constrain geophysical inversions by geology information, level-set inversion approaches

<footnote n="1">¹Laboratoire GET, Université Toulouse 3 Paul Sabatier, IRD, CNRS UMR 5563, Observatoire Midi-Pyrénées, 14 avenue Edouard Belin, F-31400 Toulouse, France. E-mail: roland.martin@get.omp.eu</footnote>

<footnote n="2">²Centre for Exploration Targeting (School of Earth Sciences), University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia. E-mail: jeremie.giraud@uwa.edu.au</footnote>

<footnote n="3">³Mineral Exploration Cooperative Research Centre, School of Earth Science, University of Western Australia, 35 Stirling Highway, Crawley, WA 6009, Australia</footnote>

<footnote n="4">⁴Laboratoire Géoressources, Equipe RING, Université de Lorraine, Rue du doyen Marcel Roubault, 54500 Vandoeuvre-lès-Nancy, France</footnote>
can be used to deform the geological structures in order to fit the geophysical data while assuming constant properties in each pre-defined geological rock unit (Zheglova et al. 2013; Li et al. 2014, 2017, 2020; Sun & Li 2016; Kadu et al. 2017; Liu et al. 2018; Zheglova et al. 2018). The number of lithologies was generally limited to three at the most, but recent works of Giraud et al. (2021a) have been able to consider an arbitrary number of geological bodies. However, as for other level-set techniques it is difficult to obtain models that are compatible with the geology. To overcome such geological drawbacks in the inversions, Giraud et al. (2022) have proposed to take into account the geology at each iteration of geophysical data inversion process by applying automated implicit geological modelling techniques to correct the inverted model and provide the most geologically compatible model as possible. This way, this geologically compatible model allows to regularize the geophysical inversion and produce a realistic and geologically acceptable model. In an on-going study recently submitted for publication (Giraud et al. 2024), this idea is further extended by integrating both the geological and geophysical information into an iterative inversion procedure in which the automated implicit geological modelling guides directly the inversion search direction. However, the physical properties remain homogeneous in each deformed rock unit. In our work, we prefer not to deform the lithologies and we rather explore another approach in which we invert the gravity data to obtain density distributions varying around the initial rock unit density values. This is done in the context of cooperative sequential inversion approach described hereafter. This allows to respect both geophysical and geological structures and allow the unit boundaries to move smoothly between adjacent lithologies and also create naturally new bodies.

To reduce uncertainties in the inverted models, two different joint geophysical inversions combining different data sets and physics are mainly used: simultaneous or cooperative inversions. As summarized in the reviews of Lelievre & Farquharson (2016), Moorkamp et al. (2016) and Ren & Kalscheuer (2020), single or multiple properties can be improved using such multiphysics and multiple geophysical data inversions approaches. This has been done for instance by inverting both gravity and gradiometric data to obtain density alone or seismic and gravity data to retrieve seismic velocities and density. Joint inversions can be performed by involving structural constraints such as structural similarities between different models as in Haber & Oldenburg (1997) and Gallardo & Meju (2003), or petrophysical constraints using clustering approaches to restrict inverted models in some space solution regions that respect some statistics about the distribution of petrophysical properties (Paasche & Tronicke 2007; Lelievre et al. 2012; Carter-McAuslan et al. 2015; Sun & Li 2015, 2016; Giraud et al. 2017, 2019a, b; Maag & Li 2018; Astic & Oldenburg 2019). Those two kinds of joint inversions allow to take advantage of the complementarity between different data sets (geophysical data and petrophysical sample data sets), but the simultaneous inversions can become computationally very expensive in CPU time and memory storage when compared to sequential and cooperative separate inversions. In this study, we thus prefer to perform such cooperative inversions, following a similar idea suggested by Lines et al. (1988) and Lines (1999) who suggested that the information and models provided by the inversion of one particular data set can be introduced as inputs and prior models for the inversion of another data set. Another reason of using cooperative inversion sequentially is that we can integrate the inversion of two different physics and their related data sets and models by running numerical tools independently: for instance seismic data and velocity and density models in one hand, and gravity data and density models in the other hand.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Integrated inversion procedure workflow: a first Taguchi-based variance analysis is applied to the prior model perturbations and an estimate of bound intervals is made to define a reduced model search space. A first ADMM-based inversion constrained by those bound intervals is performed. For each rock unit a new distribution of densities is obtained and new intervals are obtained. A second ADMM-based inversion cycle constrained by these new intervals is then performed and a final density model is obtained.</figcaption></figure>

This sequential inversion methodology is motivated by the fact that very few inversions have been applied to invert different data sets and models related to different physics (such as seismics and gravity) to reconstruct densities and seismic velocities at the crustal scale and more particularly subduction zones such as intraplate collisional areas. As a case study, we take the Pyrenees range (North Spain and South France border area) which raises economic interest due to the possible decarbonated resources potential (such as native hydrogen, Lefeuvre et al. 2022). Teleseismic first arrival times or full waveform inversion techniques have been applied to invert temporary or dense seismic arrays through different projects such as ECORS (Choukroune 1989), PYROPE (Chevrot et al. 2014, 2015) and OROGEN (Chevrot et al. 2018), and different models of the crust and the upper mantle have been provided. Gravity data inversions have been applied in that area with seismic constraints at low resolution (10 km, Dufréchou et al. 2018) at the whole scale of the Pyrenees or using density-$P$ velocity cross-gradient similarity constraints (Martin et al. 2021) at 2 km resolution along north-south dense array seismic PYROPE profile crossing the western Mauléon strong gravity anomaly. But those similarity constraints are generally used as a conservative alternative in the absence of more complete geological or geophysical prior information to derive constraints for joint inversion. This is the reason why we perform here our gravity inversions constrained by density intervals that are geophysically coherent and compatible with existing geological constraints.

In this study, we choose to perform an inversion of gravity data guided by a geologically constrained model of the Pyrenees which is available at Bureau de Recherches Géologiques et Minières (BRGM) website and is derived from the work of Wehr (2017) and Wehr et al. (2018). This model incorporates all sorts of information in the Geomodeller software, including seismic data and models. This model (also called model II and proposed by Wehr (2017, 2018) as the most representative among other models), was derived using seismic modelling for structures at depth, geological maps and also gravity forward modelling and geological knowledge. In other words, our gravity data inversions are constrained with densities restricted within regions of the solution space that respect
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Real CBA data from BGI database, (b) inverted CBA, (c) CBA computed for the prior model of Wehr *et al.* (2018) with a mantle density value of 3259 kg m⁻³ in both the Iberian and Eurasian regions or (d) for the same prior model of Wehr *et al.* (2018) but with a density mantle of 3150 kg m⁻³ in the Iberian region. Five south-north PYROPE and OROGEN seismic profiles (A–A' to E–E') and five west-east profiles (H1 to H5) superimposed on the real CBA data from BGI database. The acronyms L and SG within the figures represent Labourd and Saint-Gaudens towns locations in the Mauléon and Saint-Gaudens basins, respectively.</figcaption></figure>

geological information. These constrained gravity data inversions are a cooperative modelling effort where gravity data are used to derive intraunit density variations and to investigate more finely the hypotheses that were made in deriving the original model. A scheme illustrating the whole integrated inversion procedure is shown in Fig. 1. Our inversion workflow consists in three successive steps. In a first step, we perform a gravity data sensitivity analysis applied to perturbations of the prior model. Those perturbed prior density models are respecting the most they can the geological structures in the crust and the upper mantle. For this, we project on our discretization mesh these initial homogeneous density models related to the lithological units coming from geological information and the teleseismic data inversions, and we apply to these lithological densities a probabilistic Taguchi-based procedure (Taguchi 1987) to estimate the impact of the variations of each density unit on gravity data. This allows first estimates of density bound intervals to improve gravity data inversion. In a second step, we perform an

inversion of gravity data constrained by such intervals to obtain a first inverted density model and assess new uncertainty/bound intervals around the inverted density values. These new interval bounds are thus taken as the standard deviation (STD) distributions per lithology around the mean density value in each lithology. In a third step, gravity data are inverted constrained by these new bound intervals.

To perform such gravity inversions constrained by spatially varying density bounding constraints defined by geological units and petrophysical information, we extend and apply a specific version of the DIBC-ADMM (Disjoint Interval Bounding Constraints-based Alternating Direction Method of Multipliers) procedure of Ogarko *et al.* (2021) at each inversion step. The DIBC-ADMM is based on the ADMM methodology of Boyd *et al.* (2011) and is able to provide numerically stable inversions for inequality or bounding constraints (different and disjoint bounding intervals being possibly considered on each density voxel). Even if the DIBC-ADMM
<figure><img src="image_4.png" /><figcaption>**Figure 3.** North–south vertical sections (top, a) and west–east sections (bottom, b) of original lithologies along the seismic profiles. Units of both abscissa and vertical coordinates are in kilometres.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Horizontal sections at depth of original lithologies along the seismic profiles. Units of both abscissa and vertical coordinates are in kilometres.</figcaption></figure>

approach can be applied here, we will rather call the method as IBC-ADMM from now on everywhere in this study because we apply here an DIBC-ADMM version which introduces single Interval Bound Constraints (IBC) on each density model voxel during the inversion procedure. Our idea is to integrate bounding constraints that are first estimated through probabilistic metrics derived from the Taguchi method (Taguchi 1987; Taguchi et al. 2005) and to apply this method to pre-existing geological units.

Here, we thus aim using the IBC-ADMM bounding constraints to take into account prior geological models as well as a gradient smoothing to reduce high-frequency noisy model perturbations while allowing pre-existing geological interfaces. Besides, by taking different density anomaly intervals, it is possible to retrieve new density contrasts that can depart from the original lithological density values and expect this way to reveal the presence of hidden lithological features. But first, it is important to provide a first interval set of bounding values. This is generally enabled by classical Monte-Carlo (MC) based random probabilistic approaches. However, this may turn out to be too computationally expensive.

<table><caption>**Table 1.** Lithologies of the Pyrenees structure.</caption><thead><tr><td>Rock unit (Parameter)</td><td>Density (kg m<sup>-3</sup>)</td></tr></thead><tbody><tr><td>1. Basement/Axial zone</td><td>2660</td></tr><tr><td>2. Mesozoic sediments</td><td>2600</td></tr><tr><td>3. Cenozoic sediments</td><td>2400</td></tr><tr><td>4. Lower crust</td><td>2810</td></tr><tr><td>5. Iberian Lower Middle crust</td><td>2780</td></tr><tr><td>6. Middle crust</td><td>2755</td></tr><tr><td>7. Eurasian upper crust</td><td>2610</td></tr><tr><td>8. Hydrated mantle</td><td>3090</td></tr><tr><td>9. Upper Mantle</td><td>3260</td></tr></tbody></table>

Therefore, to mitigate this, we use a Taguchi analysis approach (Taguchi 1987) that is able to give the impact of the variation of the inversion procedure's hyperparameters represented by the density model bounds. It reduces drastically the number of runs when compared to an MC analysis as we detail in this study. Once those
<table><caption>**Table 2.** Per cent contribution $P_i$ of each lithology $i$ on the gravity data residuals for 3, 5 and 7 per cent density perturbations.</caption><thead><tr><th>Rock unit (parameter)</th><th>Minimum/maximum impact P<sub>i</sub> (per cent)</th><th>Average impact P<sub>i</sub> (per cent)</th></tr></thead><tbody><tr><td>1. Basement/Axial Zone</td><td>0.0009/20</td><td>2.4/2.46/2.46</td></tr><tr><td>2. Mesozoic sediments</td><td>0.00026/8.2</td><td>0.64/0.59/0.69</td></tr><tr><td>3. Cenozoic sediments</td><td>0.00049/1.8</td><td>0.21/0.25/0.23</td></tr><tr><td>4. Lower crust</td><td>5/10.5</td><td>7.1/7.8/7.5</td></tr><tr><td>5. Iberian Lower Middle crust</td><td>0.0012/10.1</td><td>1.6/2.6/2.7</td></tr><tr><td>6. Middle crust</td><td>2.018/20.5</td><td>5.1/4.96/5.2</td></tr><tr><td>7. Eurasian upper crust</td><td>0.00044/12.4</td><td>1.71/1.8/1.72</td></tr><tr><td>8. Hydrated mantle</td><td>8 10<sup>-6</sup>/1.8</td><td>0.1</td></tr><tr><td>9. Upper mantle</td><td>69.3/88</td><td>79</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Minimum and maximum percentage contribution $P$ (say impact for conciseness) of each lithological unit for 3, 5 and 7 per cent density variations. Curves are shown in semi-logarithmic scale to better show the differences between them.</figcaption></figure>

intervals are defined then the IBC-ADMM based inversion (see Ogarko et al. 2021, for details about the IBC-ADMM algorithm) is used to recover new models as well as new model bounding constraints consistent geologically and geophysically. The IBC-ADMM is implemented in the Tomofast-x inversion code as in Giraud et al. (2021b).

This study is organized as follows. In a first section, we introduce the gravity data and the pre-existing models that we use for the gravity data inversions at the whole scale of the Pyrenees range in south France taken as a case study to test our inversion workflow. Then, we explain the main steps of the Taguchi method and the different metrics involved without too exhaustive details and show how the density variations of the different lithologies are impacting the misfit between observed and computed gravity data. In a third section, we show the main results of our gravity inversions. In a fourth section, we discuss the models obtained and provide some mechanical and geological interpretations. More particularly we show how the method is able to recover the main features of the Pyrenean range such as the strong asymmetry at depth between the northern and southern sides of the Pyrenees as well as evidence of previously unrecognized shallow features in the Axial Zone (AZ)

and in the intraplate collision area. In the last section, we provide some conclusions and develop perspectives opened by this work.

# 2 GRAVITY DATA AND PRIOR DENSITY MODELS IN THE PYRENEES

## 2.1 Initial model of the Pyrenees

We consider the Pyrenees mountain range for validation and testing of our inversion method since it is a well-studied pilot site for developing innovative geophysical exploration methods and because some geological questions still remain (Chevrot et al. 2018, 2022). In particular, the Pyrenees region was the focus of many geological and geophysical studies because it has been the theatre of an intraplate collision defined by the subduction of the Iberian plate beneath the Eurasian Plate as described in Choukroune (1989). Furthermore, it shows a clear asymmetric structure (Chevrot et al. 2018). In the north-western side, high compressional velocities are appearing at depth with an uplift of the Eurasian upper mantle. Strong gravity anomalies appear on the Complete Bouguer anomaly (CBA) map in the Mauleon basin and Saint Gaudens
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Impact (in per cent) of 5 per cent density variations of the 9 rock units on gravity data misfit.</figcaption></figure>

[see Fig. 2 showing the gravity data provided by the BGI (Bureau Gravimétrique International), Geodesist's-Handbook 2012; Balmino et al. 2012, https://bgiobs-mip.fr]. Many studies have been led on the anomalies of that area through geological studies or geological modelling constrained by gravity data (Clariana et al. 2022) as well as inversion of seismic and gravity data inversions (Wehr 2017; Pedrera et al. 2017; Chevrot et al. 2018; Wehr et al. 2018; Pedrera et al. 2018; Lescoutre 2019; Gómez-Romeu et al. 2019; Martin et al. 2021). Furthermore seismic imaging through the inversion of teleseismic data collected by the PYROPE, IB-ERARRAY and OROGEN campaigns (Chevrot et al. 2014, 2015; Wang et al. 2016; Chevrot et al. 2018) have provided seismic velocity images (and also sections of migrated receiver functions) showing upper-mantle exhumation with an uplift close to the surface in the western and central parts of the Pyrenees. This exhumed mantle is also what we call here a Hydrated mantle. However, the reconstruction of vertical steep dipping structures is difficult when using seismic techniques. In Wehr et al. (2018), nine main lithological structures have been identified (see Figs 3 and 4), but the gravity data misfit computed with those seismic models and the geological information collected in the area still remains greater than the measurement errors (1.5 mGal approximately). This can maybe due to the presence of steep or vertical structures not easy to image

with seismic approaches or due to the presence of smaller scale bodies not taken into account in the modelling or of smaller variations in the geometry of the recovered units. As a complementary alternative to seismic imaging, we can use gravity data inversion to provide good lateral resolution of the density anomalies. We thus use the geological model consisting of 9 main geological units and the related densities obtained by Wehr (2017) and Wehr et al. (2018) at 2 km resolution as a first *a priori* model for our gravity data inversions over all the Pyrenees range. However, we want to note that the mantle densities provided by Wehr et al. (2018) seem to be overestimated in the Iberian region south from the AZ since their computed gravity anomalies are totally overestimated there. On the contrary, the Eurasian Mantle densities seem to be underestimated in the French (Eurasian) area north to the AZ. Those are also some of the main reasons leading us to new inversions in order to provide more accurate density distributions in the Iberian region and also in the AZ.

In a first step, a sensitivity analysis based on Taguchi method is applied to the prior densities of each unit for estimating the impact of variations of the unit densities on the gravity data misfit (density variations smaller than 7 per cent were considered to be enough in our analysis). This way, first estimates of reasonably physical bounding intervals are tested in our IBC-ADMM-based gravity
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Deviations of the gravity responses of all the Taguchi model combinations from the real gravity data at all observation points, for 3 per cent (left, a), 5 per cent (right, a) and 7 per cent (bottom, a) density perturbations. (b) Gravity residual Variances related to the real gravity data and the responses of all Taguchi model combinations for 3 per cent (left, b), 5 per cent (right, b) and 7 per cent (bottom, b) density perturbations. Variances are too high for the 7 per cent perturbation case.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** For a 3 per cent prior model ADMM bounding perturbation: evolution of each cost function contribution (data misfit, model and gradient smoothing, ADMM bounding constraint) to the global cost function through inversion iterations for different model damping, model gradient damping and initial ADMM parameters.</figcaption></figure>

inversions. Since shallowest densities of Mesozoic and Cenozoic are more constrained due to petrophysical information and have also smaller values than those of deeper structures (middle crust, lower crust and mantle), the interval bounds are narrower close to the surface and broader when getting deeper structures. The densities of middle/lower crust and mantle regions will thus evolve with more freedom during the inversion. As detailed in the following sections, we perform our inversions at 2 km resolution.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** For a 5 per cent prior model ADMM bounding perturbation: evolution of each cost function contribution (data misfit, model and gradient smoothing, ADMM bounding constraint) to the global cost function through inversion iterations for different model damping, model gradient damping and initial ADMM parameters.</figcaption></figure>

## 2.2 Gravity data sets

We aim at inverting the complete Bouguer gravity anomaly data provided by the BGI over a computational domain depicted in Fig. 2(a) with [-2.25,2.75]×[41.5,44.2] dimensions and covering both the French and Spain regions. This area has a dimension of 400 km in longitude × 240 km in latitude. To reduce the computational time, we make the choice here to project all those gravity data on a regular 24120 (201 × 120) data points grid that contains all the
<figure><img src="image_11.png" /><figcaption>**Figure 10.** For a 5 per cent prior model ADMM bounding perturbation and an initial ADMM weight $\alpha_{admm} = 10^{-3}$: evolution of each cost function contribution (data misfit, model gradient smoothing, ADMM bounding constraint) to the global cost function through inversion iterations for different model damping, model gradient damping and initial ADMM parameters.</figcaption></figure>

seismic PYROPE (1 and 2, i.e. A–A' and B–B') and OROGEN profiles (profiles 3, 4 and 5, i.e. C–C', D–D' and E–E', see Fig. 2). The data coordinates are converted to UTM coordinates. Average spacings around 2 km in longitude and 2 km in latitude are considered. Gravity anomalies varying between –112-and 30 mGal are observed, with the strongest anomalies being present in the western part of the Pyrenees. More particularly, the highest anomalies are clearly significant in Mauléon and Saint-Gaudens regions.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** For a 5 per cent prior model ADMM bounding perturbation and an initial ADMM weight $\alpha_{admm} = 10^{-7}$: evolution of each cost function contribution (data misfit, model gradient smoothing, ADMM bounding constraint) to the global cost function through inversion iterations for different model damping, model gradient damping and initial ADMM parameters.</figcaption></figure>

We now consider a first-order 3-D prior density model available from BRGM website at 2 km resolution and obtained by Wehr (2017) with the Geomodeller tool (Lajaunie et al. 1997; Calcagno et al. 2008; Guillen et al. 2008) which is based on a stochastic inversion procedure approach (https://www.brgm.fr/fr/logiciel/geomodeller). Geomodeller has introduced different kind of information (geological and gravity data as well as 1-D seismic profiles along six transects of PYROPE and OROGEN data sets). A first inverted density model is obtained by using Geomodeller solver and
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Histograms of inverted densities per lithology (distributions of density versus corresponding number of density values) with a ± 5 per cent variation of the prior model and final STD shown by the dashed lines.</figcaption></figure>

constrained by geological information. In a second step, a prior seismic model is obtained from this density model by using a density-velocity correlation law of Nafe–Drake type (Ludwig *et al.* 1970). Then seismic traveltimes are inverted by using a ray-tracing inversion solver detailed in Noble *et al.* (2014). Densities can be corrected again by applying the Nafe–Drake correlation law, and a new model is provided. This model differs slightly from the one published in Wehr *et al.* (2018). This can explain different means and discrepancies of 10s of mGal between computed and real gravity data, mainly in the southern/Iberian region. The different rock units and their densities are given in Table 1. This model has been obtained on a 2 km regular grid of 34800 gridpoints (240 in longitude × 160 in latitude × 36 in depth) in the Lambert93 system of coordinates. These coordinates are also converted on a UTM coordinates grid to be consistent with the gravity data UTM coordinates. The domain covered by the gravity data (400 km in E-W direction × 240 km in N-S direction) is smaller than the density anomaly model (480 km in E-W direction × 320 km in N-S direction × 72 km in depth) which is laterally extended by model invariance with an external padding to reduce much more the edge effects. These external padding regions involve two 20 km wide extra cells in N-S direction and two 30 km wide extra cells in E-W direction, one extra cell on each vertical outer layer of the density model. The whole density model has thus a total size of (480 + 2 × 30) km × (320 + 2 × 20) km × 72 km = 540 km × 360 km × 72 km. We define the density anomalies $\Delta\rho = \rho - \rho_c$ in the crust (from the upper surface down to a 30 km reference Moho depth) with a reference crust density of 2670 kg m⁻³, and $\Delta\rho = \rho - \rho_m$ beneath this Moho for a reference mantle density $\rho_m$ of 3260 kg m⁻³. In all the present study, the computations of the gravity fields involved in the forward and also

the inverse problems are performed by using a solid body approach of Blakely (1995) based on rectangular prisms as detailed in some of our recent works (Martin *et al.* 2013, 2017, 2021). We take into account the topography-bathymetry model ETOPO1 (1mn arc ≈ 1.8 km resolution) of Amante & Eakins (2009) that we interpolate on our 2 km resolution data grid. When computing the gravity data anomalies, all the rectangular prisms have the same vertical size except at the top of the computational domain where the height of each cell is adapted from the sea level up to the topography. In Fig. 2, we can observe the gravity data calculated with the *a priori* model of Wehr *et al.* (2018) on the topography (Fig. 2c) and also its difference with the observed gravity data (Fig. 2a). The strongest positive anomalies in the Pyrenean axis are located in the Labour/Mauléon (label L) region with up to 30 mGal values while strong negative anomalies up to 80 mGal are present in the eastern axial part of the Pyrenees. However, the differences between observed and prior anomalies are essentially strong in the southern and also the northern side. Indeed, in Fig. 2(c), a constant prior density of 3259 kg m⁻³ in both the Iberian and Eurasian mantle lithologies does not seem to be adequate to explain the real gravity data, more particularly in the Iberian region. Alternatively, in Fig. 2(d), the Bouguer gravity anomaly data are computed with an Iberian mantle density value of 3150 kg m⁻³ and their magnitude are closer to those of the observed anomaly data. It thus seems that a lower mantle density must be taken into account when modelling this region or in future studies in the Pyrenees. This is another reason of performing our present study. Therefore, we aim here at reducing these differences by estimating the perturbations of the prior densities that influence the data residuals the most and that will constrain the gravity data inversions. We thus perform a first
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Histograms of inverted densities per lithology (distributions of density versus corresponding number of density values) obtained after the second step (distribution between dashed lines) for interval variations of the prior models equal to ± STD (a) or ± 1.5 STD (b) for the STD obtained in the first step (same histograms as in Fig. 12). Final STD are shown by the dashed lines.</figcaption></figure>

statistical study of the impact of the variations of the prior densities to have a first estimate of the intervals of possible density values and to be able to perform inversions submitted to IBC-ADMM constraints as discussed in the following sections. Furthermore, once those intervals are defined, the IBC-ADMM method is also useful to recover intraunit variations, which is an advantage over other methods such as level set-based techniques for instance. Our geologically constrained ADMM inversion method is able to detect and estimate the density variations in both the French and Spanish parts in order to retrieve the observed anomalies as in Fig. 2(b). To show the robustness of our ADMM inversion to the bounding constraints and their consistency with the Taguchi sensitivity analysis, we perform two different inversions in Section 3: a

first one with a 3259 kg m⁻³ prior density in both the Iberian and Eurasian upper-mantle lithologies and a second one with prior density values equal to 3150 kg m⁻³ in the Iberian upper mantle (IUM) and 3259 kg m⁻³ in the Eurasian upper mantle.

## 3 DENSITY MODEL BOUNDING CONSTRAINTS ESTIMATES BASED ON TAGUCHI ANALYSIS

We aim now at improving the density model by inverting the data anomalies under bounding constraints, those constraints allowing
<table><caption>**Table 3.** Average density and STD of densities per lithology after inversion.</caption><thead><tr><th>Rock unit</th><th>Average/Std deviation with new model<br>(kg m<sup>-3</sup>)</th><th>Min/Max inverted densities<br>(kg m<sup>-3</sup>)</th></tr></thead><tbody><tr><td>1. Basement /Axial Zone</td><td>2657/64</td><td>2527/2793</td></tr><tr><td>2. Mesozoic Sediments</td><td>2593/38.4</td><td>2470/2730</td></tr><tr><td>3. Cenozoic Sediments</td><td>2405/37.4</td><td>2280/2721</td></tr><tr><td>4. Lower crust</td><td>2798/26.3</td><td>2670/2950</td></tr><tr><td>5. Iberian Lower Middle crust</td><td>2758/30</td><td>2641/2918</td></tr><tr><td>6. Middle crust</td><td>2747/37</td><td>2618/2892</td></tr><tr><td>7. Eurasian upper crust</td><td>2608/28.7</td><td>2480/2740</td></tr><tr><td>8. Hydrated mantle</td><td>3111/53.9</td><td>2935/3243</td></tr><tr><td>9. Upper Mantle</td><td>3233/24.7</td><td>3087/3411</td></tr></tbody></table>

to reduce the parameter search space. For this, we follow the global workflow of the integrated inversion procedure (Fig. 1) in which we proceed in three successive steps: -(1) a sensitivity analysis to estimate the density bounding interval constraints, followed by (2) a data inversion performed under those constraints and (3) a second inversion constrained by the intervals derived from the previous inversion. In a first step, we thus need to estimate inversion's hyperparameters defined by preliminary possible density values intervals. For this, we choose to use a statistical analysis method based on Taguchi method (Taguchi 1987; Taguchi et al. 2005). This method allows evaluating the impact of the density perturbations of the *a priori* geological units on the gravity data residuals. Generally, uncertainties around the *a priori* model can be done through a classical MC analysis. But this is very time consuming due to a too large number of random simulations (thousands of runs) to be made to explore the whole parameter space. We want here to emphasize and stress again that we use the Taguchi method as an alternative to reduce the parameter space (Mistree et al. 1993; Plazolles et al. 2015), the Taguchi method being less time consuming than MC analysis by several orders of magnitude.

The MC analysis can determine the set of parameters that are giving good results (characterized in our case by the minimized data misfit at each observation point). However, it requires a huge number of model responses that are computed for randomly generated values of input parameters. Alternatively, the Taguchi method rather computes a reduced number of computations obtained by selecting representative input parameter values and by combining them according to the so-called orthogonal arrays (OA, Taguchi 1987; R.C. Bose 1952). For instance, 9 independent computations can be performed for less than 4 input model parameters, or 27 independent computations for less than 13 parameters, etc. As a consequence, the Taguchi method determines which parameters are impacting the most the variations of the results or not. This method is thus a good candidate to estimate the uncertainties of the model densities that are expressed as probabilities of impact on the results defined here as the gravity data residuals.

In practice, the Taguchi method consists in finding how model results respond to variations of input parameters. It takes into account uncertainties around representative values of the model parameters and permits one to identify which parameter(s) have the most influence on the variation of the results of a model. Each simulation experiment of the Taguchi analysis is built by combining the representative values of the input parameters according to levels defined by OA (matrices) (Taguchi 1987; R.C. Bose 1952). This ensures that each input parameter is varied at different levels, and the effects of the parameters on the output response are studied using these OA. An OA is defined essentially by three numbers: (1) the number of levels (the number of different values taken by one model

parameter), (2) the factor number (i.e the number of columns, one column representing one parameter that can take one among the three possible values for a given experiment), and (3) the number of rows (one row per experiment, with one combination of perturbed input parameters). Here, we take a number of levels equal to 3 (perturbed input parameters), a factor number equal to 9 representative input parameters (one per lithology), and a number of rows equal to 27 (each row corresponding to one experiment defined by forward modelling runs for a combination of the 9 parameters). This number of 27 rows comes from standard OAs predefined and tabulated by Taguchi (1987, see also Phadke 2021) and is thus imposed by construction.

We want to note that the Taguchi-based variance analysis is specific to each experiment, each data set and each parameter configuration (i.e number of elements and spatial distribution of the data set, computational domain resolution/cell sizes, number of factors, 1 factor representing one rock unit and number of levels/perturbations per rock unit). For a different configuration, we will have different results. In the present study, we choose three levels (mean value and ± a perturbation) because this is what Taguchi (1987) recommends when a parameter uncertainty exists (following a Gaussian distribution or not).

In our case, the OA is the $L_{27-3}$ (i.e. 27 rows and 3 levels for each parameter). One of the main important points is to define three uncertain values for each parameter, those values corresponding to the three levels. Commonly, a normal distribution is assumed around each input parameter, with a mean value $\bar{\rho}$ and a tolerance perturbation $\Delta_{\rho}$ possibly defined by an STD. Therefore, as usually done, we choose three different representative values of each parameter, represented here by the prior density $\rho$ assigned to each of the nine lithological units, and equal to $\rho$, $\rho \pm \Delta_{\rho}$. Each representative level corresponds to one of these values (Phadke 2021). For instance : level 1 corresponds to $\rho - \Delta_{\rho}$, level 2 to $\rho$ and level 3 to $\rho + \Delta_{\rho}$.

Now, we define the following quantities that are needed for the Taguchi-based variance analysis:

(i) the number of experiments $N$ (rows of the OA $L$) here equal to 27 runs;

(ii) the number of parameters $M$ of the model (columns of the OA $L$) here equal to 9 units;

(iii) the number of degrees of freedom $f_i = L_i - 1$ is here equal to 2 for each parameter $i$ ($L_i=3$ being here the number of possible values taken by the parameter $i$);

(iv) the correction factor $CF = \frac{(\sum_{r=1}^{N} y_r)^2}{N}$, where $y_r$ is the misfit gravity data for each run $r$;

(v) the $SS_T = \sum_{r=1}^{N} y_r^2 - CF$;

(vi) the $SS_i = \frac{\sum_{j=1}^{k_i} (A_j^i)^2}{n_i} - CF$, where $k_i$ is the number of different values taken by the parameter $i$ (here $k_i$ is equal to 3), $n_i$ is the
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Maps of inverted gravity residuals' variances related to each STD-based ADMM intervals.</figcaption></figure>

<table><caption>**Table 4.** Minimum, maximum, mean STD values of gravity data misfits as well as kurtosis and skewness (no units because non-dimensional numbers) of those STD values for the inversions performed with prior density models obtained from the first inversion step and by taking density anomalies ADMM bounds defined as ± 0.5, 1, 1.5 and 2 times the STDs of this prior model.</caption><thead><tr><th>Prior model STD (kg m<sup>-3</sup>)</th><th>Min/Max (mGal)</th><th>Mean (mGal)</th><th>Data misfit STD (mGal)</th><th>Kurtosis</th><th>Skewness</th></tr></thead><tbody><tr><td>0.5</td><td>-6.3/11</td><td></td><td>0.4</td><td>41.8</td><td>1.3</td></tr><tr><td>1</td><td>-4.4/8.1</td><td></td><td>0.4</td><td>27.6</td><td>0.8</td></tr><tr><td>1.5</td><td>-3.7/7.8</td><td></td><td>0.4</td><td>31.6</td><td>2.2</td></tr><tr><td>2</td><td>-4.6/7.4</td><td></td><td>0.4</td><td>39</td><td>2.7</td></tr></tbody></table>

number of simulations (here equal to 9) performed for a given value of the parameter *i*, $A_i^j$ is the sum of the data misfit values obtained for the *j*th value of parameter *i*;

(vii) the variance $V_i = \frac{SS_i}{f_i}$, *i* being either a parameter or the error Err;

(viii) $SS_{Err} = SS_T - \sum_{i=1}^{M} SS_i$ is a metric applied on the error Err of the method;

(ix) $f_T = N - 1$ is the total number of degrees of freedom (here equal to 26);

(x) $f_{Err} = f_T - \sum_{i=1}^{M} f_i$ is the number of degrees of freedom of the error (equal to 8 = 26 - 2×9 here);

(xi) $SS'_i = SS_i - f_i \times V_{Err}$ is the pure sum square of a parameter *i*;

(xii) the per cent contribution $P_i = 100SS'_i/SS_T$ of each parameter *i* on the gravity data misfit;

(xiii) the variance ratio $F_i = V_i/V_{Err}$ for each parameter *i*.

To determine the level of confidence of the results, the per cent contribution $P_i$ and the variance ratio $F_i$ for a given parameter *i* are needed. $P_i$ measures the impact of a given parameter *i* on the response when compared to the other parameters. The higher it is, the bigger influence it has. The variance ratio $F_i$ allows one to determine the confidence level associated to a really significant contribution of an input parameter *i* to the variation of a given result (represented here by the misfit data). The $F_i$ value is then compared with the values $\alpha$ provided by the F-tables for a given
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Relative differences of the inverted models obtained respectively for a 3150 kg m⁻³ (left) or a 3250 kg m⁻³ (right) prior IUM density respect to a 3250 kg m⁻³ prior IUM density.</figcaption></figure>

level of significance. If this value is smaller than a corresponding value $\alpha$ of the F-table given in Roy (1990), then this parameter is not contributing to the variation of the response. For instance, in the present study, the number of degrees of freedom is equal to 2 for each parameter (rock unit) and equal to 8 for the error. In this case, according to the F-tables, the value $\alpha$ is equal to 3.55 for a level of significance of 0.01. For all the Taguchi analyses performed hereafter in this section for different parameter perturbations, all the variance ratios $F_i$ for each parameter $i$ are greater than a minimum value as high as $2.796 \times 10^7$. Therefore all the variance ratios $F_i$ are much greater than the significant level $\alpha$, which means that the different parameters $i$ have significant impacts $P_i$ for corresponding levels of confidence up to 99 per cent (i.e. a level of significance of 1 per cent). This validates the results of the Taguchi variance analysis in our specific case.

Here, we consider the prior densities provided in Table 1 and assigned to each geological unit of the Pyrenees model as shown in Figs 3 and 4. We then perform three different Taguchi analyses, a first one with a density perturbation $\Delta_{\rho} = 3$ per cent $\rho$ of each unit of this Pyrenees model, a second one with a perturbation $\Delta_{\rho} = 5$ per cent $\rho$, and a third one with a perturbation $\Delta_{\rho} = 7$ per cent $\rho$. We then estimate the different metrics defined previously, and more particularly the per cent contribution $P_i$ of each rock unit $i$ given in Table 2 on the gravity data misfit (i.e. difference between the observed and simulated gravity data). In Fig. 5 and Table 2, we can observe that, in average, the rock units 9 (upper mantle) and 4 (lower crust) have more influence (79 and 8 per cent average contribution, respectively) on the gravity data misfit over the whole studied area while the rock units 1 (basement/AZ) and 6 (middle crust) are presenting strong influences but on more localized areas (the AZ for rock unit 1, and the northern and southern sides apart from the AZ for rock unit 6 with

maximum per cent contributions of 20 and 20.5 per cent, respectively, and average contributions of 2.5 and 5 per cent each). Then, the Eurasian upper crust (rock unit 7) has an influence essentially in the northern part (10 per cent in average with maxima up to 12.4 per cent), while the lower middle crust (rock unit 5) has an influence in the southern/Iberian part (with a similar 10 per cent impact in average). Finally, the rock units 2 (Mesozoic sedimentary layer), 3 (Cenozoic sedimentary layer) and 8 (Hydrated upper mantle) have very few influence on the per cent contribution with less than 1 per cent average impact, and more precisely with 0.7, 0.25 and 0.2 per cent average contribution on the gravity data misfit, respectively. Those lithologies (Cenozoic and Mesozoic sedimentary basins and Hydrated upper mantle) have small impacts because the density values of the prior model are already close to the actual values and/or their volumes are significantly small compared to those of the upper mantle, the different crustal layers and the AZ. However, some maxima up to 8.2, 1.8 and 1.8 per cent, respectively, can be reached in very localized regions (see Figs 5 and 6) due to larger volumes respect to their neighbouring structures or possible greater density variations locally. On the contrary, this means also that density perturbations of all units except units 2, 3 and 8 have more influence in the inversion and then can be more prone to strong variations during inversions, and thus need to be more properly defined by inversion. Besides, the strong impact of the upper mantle on data can be explained in part because the total volume is more important than those of all the other units, this impact being still stronger outside the AZ. The upper mantle and the crustal layers have a much larger coverage in longitude and latitude, while the sedimentary Cenozoic and Mesozoic basins and the Hydrated (exhumed) upper mantle have much smaller volumes and are more localized. The AZ corresponds to a wide band at the surface with densities stronger than the sediments but its greater impact on data compared to the other
<figure><img src="image_17.png" /><figcaption>**Figure 16.** N-S sections of the difference $\rho_{\text{MOD-IBM3150}} - \rho_{\text{MOD-IBM3250}}$ between the two inverted models obtained for a $3250 \text{ kg m}^{-3}$ (first model) or a $3150 \text{ kg m}^{-3}$ (second model) IUM prior density. Thick dotted line represents the crust/mantle boundary. Negative density values show that the first inverted model overestimates density with respect to the second inverted model. Positive density values show that the first model underestimates density with respect to the second model. We applied a conservative padding value of 20 km (in N-S profiles) at the beginning and the end of each profile (white shaded area) to avoid edge effects. Note that the E-E' profile has been discarded because of too many artefacts due to poor data gravity data and possible effects of the Gulf of Lion. A-H: Aston-Hospitalet; AZ-3S: Arize-Trois Seigneurs; LA: Labourd anomaly; NPF: North Pyrenean Fault; NFPT: North Pyrenean Frontal Thrust; NPZ: North Pyrenean Zone; SG: Saint Gaudens anomaly and SPFT: South Pyrenean Frontal Thrust.</figcaption></figure>

shallow lithologies suggests that the uncertainties on its densities are also greater. Besides, for all rock units, the curves of the average or maximum per cent contributions *P* per rock unit have the same quantitative behaviour whenever 3, 5 or 7 per cent rock unit density perturbations are considered in the Taguchi analysis. However, as can be seen in Fig. 7 for the 5 per cent perturbation case, the average deviations of all the Taguchi responses combinations from the real data are more coherent physically with the strong gravity anomalies in the AZ (Saint-Gaudens and Labourd) and have smaller values in the regions apart from the AZ when compared to the 3 and 7 per cent cases. In the 5 per cent perturbation case those deviations have also intermediate values in the eastern region when compared to the

other two perturbation cases. Furthermore, the data variances obtained for the different Taguchi model combinations are too strong for 7 per cent density perturbations but are weaker for 3 per cent perturbation. Inversely, deviations to the real data are too strong for 3 per cent density perturbations but are weaker for 7 per cent perturbations. On the other hand, the data variances obtained for 5 per cent density perturbations are lying in an intermediate range between 7 and 9 mGal while collected petrophysical samples have similar density variations ranging approximately (Wehr 2017) from $2200 \text{ up to } 2500 \text{ kg m}^{-3}$ in the upper 2 km (as in the Cenozoic sediments with an average of } 2400 \text{ kg m}^{-3}$), and from 2500 up to } $2750 \text{ kg m}^{-3}$ between 4-and 6 km depths (as in the Mesozoic
<figure><img src="image_18.png" /><figcaption>Figure 17. N-S sections of the relative difference (respect to the prior model) between both inverted models obtained respectively for a 3150 kg m⁻³ prior IUM density and a 3250 kg m⁻³ prior IUM.</figcaption></figure>

sediments or the AZ with an average of 2600 kg m⁻³). Therefore, we found reasonable to select an intermediate case with 5 per cent density perturbations which correspond to density variations between ± 120 kg m³ close to the surface (as in the Cenozoic and Mesozoic sediment basins, in agreement with the shallow density logs just mentioned before), and ± 165 kg m³ in the upper mantle. We therefore decided to use these density perturbation ranges to define the IBC-ADMM bounding intervals for data inversions. In the next sections, we perform a series of inversions with ± 3 and ± 5 per cent density variations to estimate and confirm what the best choice is.

## 4 LITHOSPHERIC DENSITY MODELLING BY GRAVITY DATA INVERSION USING IBC-ADMM BOUNDING CONSTRAINTS AND DEPTH REGULARIZATION

We now proceed to a series of inversions solved by using an optimized least-square/LSQR algorithm detailed in Paige & Saunders (1982) and described in Appendix A. We first show that the choice of perturbation intervals is crucial to satisfy the convergence of the residual data inversion and the possible interval values taken by

the models at the same time. We thus perform two kinds of inversions corresponding to two different ρ perturbations (± 3 per cent and 5 per cent) of the a priori model, those perturbations defining the intervals used by the IBC-ADMM constraints as in the disjoint DIBC-ADMM version of Ogarko et al. (2021). In our case we only consider one interval per cell instead of several disjoint intervals. For sake of clarity, we recall here that we aim at minimizing the following global cost function

$$ \chi(m, D_{\text{obs}}) = \| D_{\text{obs}} - D_{\text{cal}} \|_2^2 + \alpha^2 \| W(\Delta\rho - \Delta\rho_p) \|_2^2 + \alpha_g^2 \| \nabla(\Delta\rho) \|_2^2, \quad (1) $$

under the constraints of interval bounds as defined in Appendix A (see eq. A6). Δρ is the density anomaly and Δρₚ the prior density anomaly relative to a reference density model ρₚₐ as defined in Appendix A. Dobs and Dcal are the measured and the calculated Bouguer gravity data using rectangular prisms as in the GBOX code approach of Blakely (1995). The parameters α and αg are the damping and smoothness regularization weights applied to the density model and the model gradients, respectively. In summary, the inversion constrained by interval bound constraints consists in minimizing χ(m, Dobs) such that the density anomaly mk = Δρk of each kth model cell belongs to a set of possible intervals Bk defined
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Left: receiver functions showing similar geological interfaces (see figures from Chevrot *et al.* 2018) as in (right) the prior and inverted models at 2 km resolution along the 5 south-north seismic PYROPE vertical profiles. Black boxes are showing the subduction zones of interest. Units of both abscissa and vertical coordinates are in kilometres. A 3250 kg m⁻³ IUM prior density is considered.</figcaption></figure>

as

$$B_k = \bigcup_{l=1,L}^{L} [a_{kl}, b_{kl}] \text{ with } 1 \le k \le n \quad (2)$$

where $B_k$ is the interval binding each of the $n$ model cells. For the $k$th model cell, $B_k$ is defined as a union of disjoint intervals where $a_{kl}$ and $b_{kl}$ are the lower and upper bounds for the $k$th model cell and for each $l$th geological lithology unit. In the more general disjoint IBC-ADMM (i.e DIBC-ADMM) method, $L^k$ is the total number of possible units per cell. But in our single interval IBC-ADMM version used here, we assume that each model cell can belong to only one unit. Therefore, as $L^k=1$ for each cell $k$ here in our study, $B_k$ is expressed by only one interval per cell as

$$B_k = [a_k, b_k] \text{ with } 1 \le k \le n \quad (3)$$

These interval bound constraints are applied through the minimization of the ADMM bounding constraint cost function given by eqs (A6)-(A9) in Appendix A, to which we refer the reader to for more information.

In Figs 8 and 9, the evolution of the different cost functions are shown for different IBC-ADMM weights and without damping weights on the model. It can be observed that inversions are converging much faster with exponential decrease of the ADMM cost to very low values for the ± 5 per cent perturbation case when compared to the ± 3 per cent case. However, in all cases (see also Figs 10 and 11), only initial IBC-ADMM $\alpha_{admm}$ weight values strictly smaller than $10^{-3}$ are allowing to make all the cost functions to decrease drastically.

In summary, after running the different series of inversions, we only retain the case of model perturbation bounds of ± 5 per cent

because those perturbation bounds are the most reasonable choice to make and are corresponding to maximum ≈ ± 160 kg m⁻³ variations of the densities, which are lying in an acceptable and maximum variation range between adjacent layers.

In all inversions, the LSQR algorithm (Paige & Saunders 1982) is used (see eqs in Appendix A). 100 global outer loop LSQR inversion cycles are run with 200 inner loop iterations at each LSQR inversion cycle to reach satisfactory convergence and to obtain a normalized data misfit (respect to the norm of the real data) just below 0.01 per cent. Besides, a Haar wavelet compression of 10 per cent of the sensitivity kernel is used as in Martin *et al.* (2013) and Ogarko *et al.* (2024) to reduce the computation requirements (without compression the kernel matrix had a size of ≈ 24120 data × 1.4113 10⁶ [(Nlong = 242) × (Nlat = 162) × (Ndepth = 36)] model grid cells leading to ≈ 34.041 billion elements). Both the padding stretching function and the compression procedure allow for a significant reduction in computational time. Furthermore, since the models can still be important in size and that many runs have to be done to test the different weights and different inversion configurations, further computing accelerations are needed. Therefore, a parallel version has been designed including the compression of the kernel matrix. All the inversions are running on 400 Skylake Intel processors in less than 3 hr on a supercomputing platform located in the CALMIP mesocentre of Toulouse (France). To perform such inversions we use the parallel open source inversion code Tomofast-x (read Giraud *et al.* 2021b, 2023; Ogarko *et al.* 2024, for more details about the algorithms used and the functionalities of this code).

Due to non uniqueness of the solutions obtained by the inversions, it is crucial to constrain the minimization of the
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Inverted density models at 2 km resolution along the first four south-north seismic PYROPE vertical profiles. The different interfaces of interest are shown. Units of both abscissa and vertical coordinates are in kilometres. A 3250 kg m⁻³ IUM prior density is considered.</figcaption></figure>

global cost function given in eq. (1) by adding supplementary regularization cost functions (see also Appendix A for further details). It is thus necessary to define the reasonable range of the different regularization weighting parameters (model damping weight $\alpha$, gradient damping weight $\alpha_g$ and starting IBC-ADMM weight $\alpha_{admm}$ applied to the bounding interval constraints as in appended eq. A6) that are related to each regularization cost function. For this purpose, a reasonable number of inversions (less than 20 inversions) with different constraints weights have been tested for the 5 per cent density perturbation IBC-ADMM interval case. For each inversion run corresponding to a given set of regularization parameters ($\alpha, \alpha_g, \alpha_{admm}$), we analyse the evolution and convergence of the different cost functions contributing to the global cost function as shown in Figs 10 and 11. Besides, in order to constrain properly the models inside the IBC-ADMM bounding intervals, it is important to decrease the IBC-ADMM cost function to low values during the inversion process. To achieve this, we choose a relatively low initial IBC-ADMM weight ($10^{-8}$) and increase it by a factor 2 at those inversion iterations for which the normalized gravity data misfit term (non-dimensional) decreases under a given threshold taken here as $10^{-4}$. Subsequently, both data and ADMM costs are decreasing in a balanced way and are converging to small values as being desired to obtain bounded density models that explain the geophysical data at the same time. The same kind of inversion (not shown here for sake of clarity) was also done with an increase by a smaller factor of 1.1 to avoid too sharp variations of the IBC-ADMM constraints, but the results were very similar. Therefore, all the results shown here are obtained with an increase

of the initial IBC-ADMM weight by a factor 2, which is a reasonable increasing factor.

The models are obtained for different values of the model damping weight $\alpha$ between $10^{-8}$ and $10^{-4}$, of the gradient damping weight $\alpha_g$ between $10^{-11}$ and $10^{-7}$, and for a starting IBC-ADMM weight $\alpha_{admm} \le 10^{-7}$. For a model damping weight $\alpha > 10^{-4}$, a gradient damping weight $\alpha_g > 10^{-9}$ and a starting IBC-ADMM weight $\alpha_{admm} > > 10^{-7}$, the different terms of the global cost function are not decreasing significantly under reasonable normalized cost thresholds of 0.01 per cent and inversions are not converging. For too small values of gradient and model damping weights, solutions are not well constrained and are showing too many high-frequency numerical artefacts more specifically close to the surface. For gradient damping weights $\alpha_g$ between $10^{-10}$ and lower values, the solutions are very similar and high-frequency artefacts are still appearing close to the surface and in the crust, and structures at depth are not well constrained with too large wavelength models from the surface down to the mantle region. Therefore, we only retain the more reasonable model damping, gradient damping and starting IBC-ADMM constraints weights that give physical results without too many artefacts and that correspond to the inflection points of L-curve type as depicted in Hansen (2001). Therefore, after analysing the series of runs of Fig. 11, we can consider a model weight $\alpha = 10^{-4}$, a gradient damping weight $\alpha_g = 10^{-9}$, and a starting IBC-ADMM weight of $10^{-7}$ as optimal regularization weight values. Among all the inversion runs, the optimal models are obtained for those regularization parameters.
<figure><img src="image_21.png" /></figure>

**Figure 20.** N-S sections of the difference between the inverted and prior models for a 3250 kg m⁻³ IUM prior density (2 km resolution along the first four south-north seismic PYROPE vertical profiles). Thick dotted line represents the crust/mantle boundary. Units of both abscissa and vertical coordinates are in kilometres. Negative density values show that the prior model overestimates density with respect to the inverted model. Positive density values show that the prior model underestimates density with respect to the inverted model. We applied a conservative padding value of 20 km (in N-S profiles) at the beginning and the end of each profile (white shaded area) to avoid edge effects. Note that the E-E' profile has been discarded because of too many artefacts due to poor data gravity data and possible effects of the Gulf of Lion. A-H: Aston-Hospitalet; AZ-3S: Arize-Trois Seigneurs; LA: Labourd anomaly; NPF: North Pyrenean Fault; NFPT: North Pyrenean Frontal Thrust; NPZ: North Pyrenean Zone; SG: Saint Gaudens anomaly and SPFT: South Pyrenean Frontal Thrust.

In Figs 12 and 13 and Table 3, the distribution of the number of voxels and their related densities obtained after inversion are shown for each original lithology. We observe clear regular density changes and distribution shapes for all the original lithologies with minimum STDs of ±24.7 kg m⁻³ for unit 9 (upper mantle) and strongest STDs of ±64 kg m⁻³ for unit 1 (AZ) and ±53.9 kg m⁻³ for unit 8 (hydrated mantle). However, for units 4 (Lower crust) and 9 (upper mantle) a bi-modal distribution is observed, with two peaks of 2780 and 2810 kg m⁻³ in unit 4 and two more individualized peaks of 3210 and 3250 kg m⁻³ in unit 9 (see the STDs in Table 3). This can explain also that in the two different models I and II of the Pyrenees proposed in Wehr (2017) and Wehr et al. (2018), the upper-mantle density is 3200 kg m⁻³ in model I and 3259 kg m⁻³ in model II. This also means that, in the regions of the upper mantle and the lower crust, significant variations and increases of mass are appearing in some parts of the model. This is in agreement with the strong impacts of the density variations of these units on the gravity data (see previous section on Taguchi analysis).

In a second step, we perform a new series of inversions with four different new IBC-ADMM interval sets based on the STDs of the models obtained in the first inversion step. To see the impact of the *a priori* model on the inverted models, we rebuild a new *a priori* model defined by the inverted model obtained previously, and we define new IBC-ADMM bound intervals using the STDs (see

Table 3) obtained previously in the first step. After different tests, we have estimated that these new inversions can be performed with the same model and gradient damping weights as well as the same starting IBC-ADMM weights used in the first step.

In this second inversion, we tested four different sets of IBC-ADMM bound constraints (using 0.5, 1, 1.5 or 2 STD IBC-ADMM intervals). All the inverted models are reproducing the main interfaces of the prior model except in the Ebro basin southern side of the westernmost profile passing through the Mauleon Basin where strong and low-density anomalies are appearing there. However, it is possible to discriminate among these four inverted models those that can be kept. Indeed, as can be observed in Fig. 14, the global variance of the data misfit is computed at each data point and for the four IBC-ADMM STD constraints, and is lower than 2.5 mGal for the 1, 1.5 and 2 STD cases which have similar data residual variances. When computing the differences between the data and its mean in a window around each data point, the maximum of those differences are also lower than 2.5 mGal for these three best cases. The worst case is for 0.5 STD case which shows residual variances greater than 4 mGal in an area located south to Saint-Gaudens (corresponding to SG acronym on Figs 1 and 14). In all cases, data residual variances are approximately equal to 0.4 mGal and are of the order of the noise of the real data (≈ 2 mGal). Due to poor data quality in some areas, some data residuals are above the data noise,
<figure><img src="image_22.png" /><figcaption>**Figure 21.** Vertical sections of prior (Left) and inverted (right) models at 2 km resolution along the five west-east vertical profiles. Units of both abscissa and vertical coordinates are in kilometres. A 3250 kg m⁻³ IUM prior density is considered.</figcaption></figure>

and inverted data can not fit the real data accurately as in the 0.5 STD case with residual values reaching up to 11 mGal values.

In Table 4, the maximum of the absolute local variances reaches the lowest value for the 1 and 1.5 STD cases. Furthermore, kurtosis and skewness of all the data residual maps have the lowest values in the 1-STD case. Therefore, in the next section, we will now just consider the 1-STD case for discussion and interpretation of the inverted models.

Finally, in order to show the robustness of the IBC-ADMM algorithm, we performed a second inversion with similar regularization parameters as for the previous inversions and with a same prior model except in the IUM where ambiguities seem to be. In the IUM-, we assign a prior density of 3150 kg m⁻³. As can be observed in Figs 15 and 16, we obtain after inversion a very similar distribution of densities even in the IUM, a little increase of the densities in the AZ and in the Mesozoic and Cenozoic sedimentary units. However, the same geological interpretations and conclusions can be made and very similar results are obtained through these two inversions (as can be observed when comparing north-south and east-west profiles of Figs 20 and 23 for a 3250 kg m⁻³ prior IUM density and those shown in Figs A1 and A2 provided in Appendix A for a 3150 kg m⁻³ prior IUM density). More specifically, as can be observed in Fig. 17, very small relative differences between both inverted models are observed everywhere except in the AZ which is subject to relative differences up to 1 per cent inside but greater than 4 per cent locally due to the possible presence of denser bodies

inside. Differences close to 1 per cent can also be distinguished at the boundaries of the sedimentary basins, of the AZ and at the boundaries between the coupled IUM/lower crust system and the Eurasian upper mantle/lower crust system. These relative differences between both inverted models remain small, are geologically consistent and/or are localized at the boundaries between rock units where ambiguities are inherent by nature. Furthermore, in Fig. 15, the differences between both inverted models and the prior model II from BRGM are very similar with negative relative differences of approximately up to 2–4 per cent in the IUM and up to 1 per cent in the Iberian crust. Besides, positive relative differences up to 1–3 per cent in the AZ and up to 1 per cent in sedimentary basins, Eurasian lower crust and mantle are obtained. This shows that the variations of the rock units’ densities and their impacts on data estimated by the Taguchi method are closely and consistently related. After both inversions, it appears that both IUM and Iberian lower crust seem to extend less to the north than the prior model (called ‘mantle shortening’ on Figs 16 and 20). However, further investigation will be needed to confirm this geological and geomechanical feature.

## 5 RESULTS OF THE INVERSIONS AND INTERPRETATIONS

For sake of simplicity, we only consider now the inverted model obtained for a 3250 kg m⁻³ prior IUM density, since the two prior
<figure><img src="image_23.png" /><figcaption>**Figure 22.** Inverted models at 2 km resolution along the five east-west seismic profiles. The different interfaces of interest are shown. Units of both abscissa and vertical coordinates are in kilometres. A 3250 kg m⁻³ IUM prior density is considered.</figcaption></figure>

IUM densities (3150 or 3250 kg m⁻³) provide very similar inverted models. Figs 18–20 show five different south-north oriented vertical profiles after inversion. Figs 21–23 in turn display five west-east profiles (H1 to H5) passing along or parallel to the main axis of the Pyrenees (line H2) and covering the axial, southern and northern part of the Pyrenees. The sections from Figs 18–20 clearly reproduce the same trends as the interfaces depicted by the seismic receiver functions (Fig. 18a) obtained along the similar seismic profiles A–A' to E–E' of PYROPE and OROGEN surveys used in the studies of Chevrot *et al.* (2018). Fig. 20 also clearly evidences the slight overestimation (maximum anomalies up to approximately 100 kg m⁻³) respectively away from the AZ of the densities of unit 5 in the prior model (Iberian lower-middle crust) while a mean density of unit 7 (Eurasian plate’s upper crust) tends to be slightly underestimated (maximum positive and negative anomalies up to approximately +50 and −50 kg m⁻³, respectively). This is evidenced by the subvertical contrast between south and north at depth (below 30 km). Besides, the eastern part of the Pyrenean belt has been widely affected by Neogene opening of the Gulf of Lion. Although mechanisms are still debated, this geodynamic event caused abrupt changes in the crustal geometry, thermal state and

lithosphere–asthenosphere boundary (LAB) geometry (e.g. Fullea *et al.* 2021; Gunnell *et al.* 2008; Huyghe *et al.* 2020; Jolivet *et al.* 2020; Torné *et al.* 2023). We believe these structural changes may cause anomalous modelling results if they are not well constrained. Furthermore, the lack of density data in the most eastern part does not allow us to conduct the inversion with an acceptable level of confidence. We thus decided to discard profile E–E’ from further interpretation due to suspected edge effects and lack of data around this section that can induce non geological artefacts on the profile. Both Figs 20 and 23 are displaying well the main features of the intraplate collision of the Iberian and Eurasian plates and associated subduction processes (when present) with negative density anomalies reaching 30–40 km depth in the western and central parts of the Pyrenees, and rising of dense lower crust and upper-mantle material (strong positive density anomalies) up to 10 km depth mainly in the AZ part (H2, H3 and H4 profiles). However, towards the eastern part, where the plate subduction is not indicated by seismological data and a thinning and flat Moho between 20 and 30 km depth is present, a slightly lesser dense mantle is evidenced apart the AZ essentially in the Iberian region (H3 and H4 profiles) and also a little bit in the AZ (H2 profile). Note that the Iberian slab plunges
<figure><img src="image_24.png" /><figcaption>**Figure 23.** E-W sections of the difference between the inverted model and the prior model for a 3250 kg m⁻³ IUM prior density. Thick dotted line represents the crust/mantle boundary. Negative density values show that the prior model overestimates density with respect to the inverted model. Positive density values show that the prior model underestimates density with respect to the inverted model. We applied a conservative padding value of 40 km (in E-W profiles) at the beginning and the end of each profile (white shaded area) to avoid edge effects. Note that the E–E’ profile has been discarded because of too many artefacts due to poor data gravity data and possible effects of the Gulf of Lion. A-H: Aston-Hospitalet; AZ-3S: Arize-Trois Seigneurs; LA: Labourd anomaly; NPF: North Pyrenean Fault; NFPT: North Pyrenean Frontal Thrust; NPZ: North Pyrenean Zone; SG: Saint Gaudens anomaly and SPFT: South Pyrenean Frontal Thrust.</figcaption></figure>

at a maximum depth of 60 km (Fig. 20), which is much less than the 70 km imaged in the works of Chevrot *et al.* (2018). As can also be observed in the four north–south oriented vertical profiles (A–A’ to D–D’), the densities of the AZ close to the surface should be increased up to ± 60 kg m⁻³ when compared to the prior model. This can also be observed more particularly in the middle of the profiles passing through or just aside from the Saint-Gaudens anomaly region (B–B’ and C–C’ profiles).

Besides these first order observations, when compared to the 2 km resolution model of Martin *et al.* (2021), the results highlight smaller scale structures exhibiting negative and positive density anomalies (max ± 250 kg m⁻³) mainly in the middle and upper crust from the AZ. The BB’ profile shows for example the presence of denser material in the upper 10 km of the model and probably linked to the rooting of the Maladeta plutonic complex, also observed on the H4 latitudinal profile and a little bit on the H3 profile as can be observed on Figs 20 and 23 (Clariana *et al.* 2022). BB’ profile also

nicely displays the underestimation of the triassic evaporites linked with the Rialp thrust sheet and linked with the southverging stacking of AZ crustal units as in Clariana *et al.* (2022, see Fig. 20). Besides, on that same BB’ profile, there is also a negative density difference below the shallow dense density anomaly of Saint-Gaudens as can be observed on that same figure. This suggests that the dense body responsible of this gravity anomaly is disconnected from the underlying mantle, and is more likely related to a tectonic lense detached from the mantle into the crust (Casas *et al.* 1997; Mouthereau *et al.* 2014; Angrand *et al.* 2020, 2022). We can also suspect that the slab is shorter and plunges less deep as in the prior model. Another possibility is that the slab has a steeper slope geometry. In addition, both this profile and profile C–C’ show negative values shallower than the décollement surface in the Iberian crust. This is an interesting feature that may be explained either by less dense material, either upper crust or sediment, accreted at the footwall of the antiformal stack, or a different Iberian margin geometry than
<figure><img src="image_25.png" /><figcaption>**Figure 24.** Horizontal sections of inverted and prior models at different depths (10–50 km) and their differences at 2 km resolution.</figcaption></figure>

the one modelled on the BRGM prior model. AA' profile shows positive anomalies linked to the Labourd uplifted subjacent mantle body (Ford et al. 2022) and other positive anomalies further north in the Aquitaine basin linked to a dense body in the upper crust whose origin still remains to be assessed (Angrand et al. 2018). On CC' profile, the inverted model evidences the up to 10–15 km roots of the Aston Hospitalet metamorphic domes and further south for the Mont Louis granite on latitudinal H3 profile, and also a little bit on H4 profile (Cochelin et al. 2017). On DD' profile denser bodies are suspected at the southern edge of the AZ under the Jurassic to Cenozoic cover potentially linked to major detachment faults (Ford et al. 2022). On the latitudinal H2 profile (Fig. 23), we can see that strong density anomalies are appearing in Saint-Gaudens anomaly region from the upper surface down to 15 km approximately. In the eastern part of those profiles, the lower crust and the Moho are becoming flat without dense mantle material in the middle and lower crust. Another interesting feature is the mass excess of the inverted model respect to the prior model below the Ebro basin (profile A–A' of Fig. 20). Although too deep to be explained by sedimentary processes (down to 20 km depth) this may reveal

strong heterogeneities within the crust, due to complex geodynamic history, today masked by Ebro sedimentary basin. This is in line with recent work of Asti et al. (2022), which proposed that the Ebro basin was an active part of the Iberia-Europe transtensive system. Within the upper crust, these longitudinal profiles show mass excess below the Mouthoumet and Basques Massifs (see profiles H1 and H3 in Fig. 23). These higher density values in the inverted model compared to the prior model may reveal a too simple geometry of the North Pyrenean Massif in the BRGM model. Indeed, the Basques Massifs plunge and continue eastward below the NPZ. The lack of density below the Mouthoumet Massif may be due to crustal units rooting deeper than previously modelled.

Our method is thus able to provide more details in some areas that have been less investigated with geophysical methods.

In Figs 24, we can also see that at depths greater than 32 km (32, 40 and 50 km) the differences between prior and inverted models are significant along the axis of the Pyrenees and more particularly in the Iberian side with lower density values than in the French part. Negative variations of −100 kg m⁻³ are appearing in the Iberian side and positive variations of +50 kg m⁻³ in the AZ and the
French part. These lower density values in the Iberian lithospheric mantle may originate from mantle depletion during the Variscan Orogeny (e.g. Torné et al. (2023) and references therein). This is also in agreement with the impact percentage $P_i$ of lithological units 4 (the lower crust) and 9 (the upper mantle) on the gravity data. The lower crust and the upper mantle are thus mainly varying at depth after inversion when compared to the other lithologies (see also Table 2). The prior model densities provided by Wehr (2017) were overestimated at depth (from the lower crust down to the upper mantle) in the Spanish part and underestimated in the French part. In the Iberian side, an averaged mantle density close to $3150 \text{ kg m}^{-3}$ should be more adequate than the constant density value of $3259 \text{ kg m}^{-3}$ proposed by Wehr (2017) and Wehr et al. (2018) in both the Iberian and Eurasian mantles. The new mantle density estimates obtained with our method are more coherent with the $Vp$ and $Vs$ seismic velocities obtained by using tomography techniques (Chevrot et al. 2022), such velocities being smaller in the Iberian side than in the Eurasian (north) side or the IUM being deeper (outside the main Pyrenees axis) than the Eurasian mantle. Both the Iberian LAB depth and Iberian lithospheric mantle have indeed been shown to be highly heterogeneous below the Iberian Peninsula basins (Fullea et al. 2021; Torné et al. 2023). LAB depths range from 190 to 200 km (below the Pyrenees) to less than 70 km depth in the Western Mediterranean Neogene (Torné et al. 2023).

During inversion, new structures are thus appearing and densities vary within each initial/prior rock unit. An interesting future research avenue would be to use these variations to define new rock units for inversion and to invert for both the densities and the different units boundaries in a same inversion framework. Both the geometry of each rock unit and the density variations within each unit could thus be recovered, instead of just changing the densities within. This could be done in a more global iterative and sequential procedure by introducing implicit geological constraints at each inversion cycle. For instance, a deformation of the geometry of each rock unit could be obtained (Giraud et al. 2023, 2024) to fit both the gravity data and the geology. And this could be done by using level set inversions to investigate alternative structural scenarios as in Giraud et al. (2021a, 2023). This could then be introduced as a new prior model in the IBC-ADMM inversion to obtain heterogeneous densities within each rock unit. In the mid/long term, it would thus be possible to obtain a full model geophysically and geologically consistent.

# 6 CONCLUSION

In this study, we have developed a probabilistic Taguchi-based methodology to estimate the impact of the model density variations assigned to each rock unit on the gravity data misfit. This has allowed to identify the rock units that primarily contribute to the gravity data for a specific geological or geophysical model. We have shown that the method can also be used to obtain an estimate of the interval density model bounds for each model cell and thus to constrain the gravity data inversion by those bound constraints. This method is general in nature and is readily applicable to other linear or quasi-linear geophysical inversion techniques such as magnetic or traveltime seismic data inversion, or can also be easily extended to nonlinear inversion methods (full waveform, magneto-telluric, electrical resistivity data, etc.). It allows to address a knowledge gap by estimating uncertainty in deterministic settings. It could be applied not only to the densities of each unit as done here but also more generally to each cell of the whole computational domain. Using the Taguchi-derived bounds, we apply the IBC-ADMM to

the inversion of gravity data over the Pyrenees range under the constraint of a prior geophysical model that is consistent geologically. We have shown that the gravity data are mainly sensitive to the density perturbations in the Pyrenees AZ, in the upper mantle as well as in the lower and middle crust. More particularly, after inversion, the following significant differences in density were observed in comparison to the prior geological model:

(i) a slight increase in crustal densities is observed in the western part of the intraplate collisional zone, particularly within the first 10 km depth close to the topography;

(ii) from the middle crust down to the upper mantle, we observe a decrease of density equal to approximately $100 \text{ kg m}^{-3}$ in the southern (Iberian) side and an increase of $50 \text{ kg m}^{-3}$ in the northern (Eurasian) side as well as in the AZ close to the surface-;

(iii) a possible shortening of the lower Iberian crust subducted below the lower Eurasian crust in western AZ location.

Besides, questions are still remaining about the effects of these density variations on the topography. Indeed, it has been proposed that the Iberian and Pyrenean topography is dynamically supported by density changes in the sub-lithospheric mantle since the Neogene (e.g. Huyghe et al. 2020; Conway-Jones et al. 2019; Fullea et al. 2021). An interesting future of the present work would be to question these new findings about densities variations regarding variations of topography with local isostatic modelling. Furthermore, these results and their interpretations could pave the way to further inversions improved at finer discretization resolutions on more localized areas in the Pyrenees or other regions. Inversions could be constrained not only by better bounding interval constraints derived from petrophysical data measurements and/or accurate seismic models coming from dense seismic acquisitions, but also by geometrical constraints using level sets.

# ACKNOWLEDGMENTS

We thank the Bureau de Recherche Géologique et Minière (BRGM) who provided the models published in Wehr (2017) and Wehr et al. (2018). We thank also the supercomputing project P1138 at CALMIP mesocenter in Toulouse France for the supercomputing hours needed for this study. VO thanks MinEx for its CRC grant. This is MinEx CRC document no. 2023/33. JG acknowledges support from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 101032994. The authors express also their consideration to Sébastien Chevrot and Lucia Seoane from CNRS for interesting discussions about the modelling and history of the Pyrenees and for the gravity data set extracted from the database of the BGI (Bureau Gravimétrique International), respectively. We acknowledge Anne-Marie Cousin from CNRS for her precious expertise and help in improving the figures quality. Finally, we are grateful to Guillaume Caumon for fruitful discussions about constrained inversions.

# AUTHOR CONTRIBUTIONS

Roland Martin, Vitaliy Ogarko and Jérémie Giraud are the main contributors to this study. Bastien Plazolles provided the expertise on Taguchi analysis dispersion analysis and participated to the manuscript writing. Sonia Rousse, Paul Angrand and Melina Mancouin provided expertise on Pyrenees range geology interpretation and participated to the manuscript writing.

All authors reviewed the manuscript.
## COMPETING INTERESTS

The authors declare no competing interests.

## DATA AVAILABILITY

The Bureau de Recherche Géologique et Minière (BRGM) provided us the voxet database with rock units related to the geological models published in Wehr (2017) and Wehr et al. (2018) and downloadable from Geomodeller website on demand to the BRGM (https://www.brgm.fr/fr/logiciel/geomodeller). The BGI provided us the gravity data set used here for our inversions (https://bgiobs.mip-fr/).

## REFERENCES

*   Amante, C. & Eakins, B.W., 2009. ETOPO1 1 arc-minute global relief model: procedures, data sources and analysis., in *NOAA Technical Memorandum NESDIS, NGDC-24*, Vol. 24, p. 19, https://repository.library.noa.a.gov/view/noaa/1163.

*   Angrand, P., Ford, M. & Watts, A.B., 2018. Lateral variations in foreland flexure of a rifted continental margin: the aquitaine basin (sw france), *Tectonics*, **37**(2), 430–449.

*   Angrand, P., Mouthereau, F., Masini, E. & Asti, R., 2020. A reconstruction of iberia accounting for western tethys–north atlantic kinematics since the late-permian–triassic, *Solid Earth*, **11**(4), 1313–1332.

*   Angrand, P., Ford, M., Ducoux, M. & De Saint Blanquat, M., 2022. Extension and early orogenic inversion along the basal detachment of a hyper-extended rifted margin: an example from the central Pyrenees (France), *J. geol. Soc.*, **179**(1), jgs2020–003.

*   Asti, R., Saspiturry, N. & Angrand, P., 2022. The mesozoic iberia-eurasia diffuse plate boundary: A wide domain of distributed transtensional deformation progressively focusing along the north pyrenean zone, *Earth-Sci. Rev.*, **230**, 104040, doi:10.1016/j.earscirev.2022.104040.

*   Astic, T. & Oldenburg, D.W., 2019. A framework for petrophysically and geologically guided geophysical inversion using a dynamic Gaussian mixture model prior, *J. geophys. Int.*, **219**(3), 1989–2012.

*   Balmino, G., Vales, N., Bonvalot, S. & Briais, A., 2012. Spherical harmonic modelling to ultra high degree of Bouguer and isostatic anomalies, *Journal of Geodesy*, **86**(7), 499–520.

*   Blakely, R.J., 1995. *Potential Theory in Gravity and Magnetic Applications*. Cambridge University Press, Cambridge, UK.

*   Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J. et al., 2011. Distributed optimization and statistical learning via the alternating direction method of multipliers, *Found. Trends Mach. Learn.*, **3**(1), 1–122.

*   Calcagno, P., Chilès, J., Courrioux, G. & Guillen, A., 2008. Geological modelling from field data and geological knowledge. Part I. Modelling method coupling 3D potentiell-field interpolation and geological rules, *Phys. Earth planet. Inter.*, **171**, 147–157.

*   Carter-McAuslan, A., Lelièvre, P. & Farquharson, C., 2015. A study of fuzzy c-means coupling for joint inversion, using seismic tomography and gravity data test scenarios, *Geophysics*, **80**, w1–w15.

*   Casas, A., Kearey, P., Rivero, L. & Adam, C.R., 1997. Gravity anomaly map of the pyrenean region and a comparison of the deep geological structure of the western and eastern pyrenees, *Earth planet. Sci. Lett.*, **150**, 65–78.

*   Chevrot, S., Villaseñor, A. & Sylvander, M. & the PYROPE Team, 2014. High-resolution imaging of the Pyrenees and Massif Central from the data of the PYROPE and IBERARRAY portable array deployments, *J. geophys. Res.*, **119**(8), 6399–6420.

*   Chevrot, S., Sylvander, M., Diaz, J., Ruiz, M., Paul, A. & the PYROPE Working, group, 2015. The Pyrenean architecture as revealed by teleseismic P-to-S converted waves recorded along two dense transects, *Geophys. J. Int.*, **200**, 1096–1107.

*   Chevrot, S. et al., 2018. The non-cylindrical crustal architecture of the Pyrenees, *Sci. Rep.*, **8**, p. 9591.

*   Chevrot, S. et al., 2022. Passive imaging of collisional orogens: a review of a decade of geophysical studies in the pyrénées, *BSGF - Earth Sci. Bull.*, **193**, 1, doi:10.1051/bsgf/2021049.

*   Choukroune, P., 1989. The ecors pyrenean deep seismic profile reflection data and the overall structure of an orogenic belt, *Tectonics*, **8**(1), 23–39.

*   Clariana, P. et al., 2022. Basement and cover architecture in the central pyrenees constrained by gravity data, *Int. J. Earth Sci.*, **111**, 641–658.

*   Cochelin, B., Chardon, D., Denèle, Y., Gumiaux, C. & Le Bayon, B., 2017. Vertical strain partitioning in hot variscan crust: Syn-convergence escape of the pyrenees in the iberian-armoricán syntax, *Bull. Soc. géol. Fr.*, **188**(6), 39, doi:10.1051/bsgf/2017206.

*   Conway-Jones, B.W., Roberts, G.G., Fichtner, A. & Hoggard, M., 2019. Neogene epeirogeny of iberia, *Geochem. Geophys. Geosyst.*, **20**(2), 1138–1163.

*   Dufréchou, G., Tiberi, C., Martin, R., Bonvalot, S., Chevrot, S. & Seoane, L., 2018. Deep structure of pyrenees range (sw europe) imaged by joint inversion of gravity and teleseismic delay time, *J. geophys. Int.*, **214**(1), 282–301.

*   Ford, M. et al., 2022. Evolution of a low convergence collisional orogen: a review of pyrenean orogenesis, *BSGF - Earth Sci. Bull.*, **193**, 19, doi:10.1051/bsgf/2022018.

*   Fullea, J., Negredo, A., Charco, M., Palomeras, I., Afonso, J. & Villaseñor, A., 2021. The topography of the iberian peninsula from integrated geophysical-petrological multi-data inversion, *Phys. Earth planet. Inter.*, **314**, 106691, doi:10.1016/j.pepi.2021.106691.

*   Gallardo, L.A. & Mejú, M.A., 2003. Characterization of heterogeneous near-surface materials by joint 2d inversion of dc resistivity and seismic data, *Geophys. Res. Lett.*, **30**(13), doi:10.1029/2003GL017370.

*   Geodesist's-Handbook, 2012. The international gravimetric bureau, *J. Geod.*, **86**(10), 946–949.

*   Giraud, J., Pakyuz-Charrier, E., Jessell, M., Lindsay, M., Martin, R. & Ogarko, V., 2017. Uncertainty reduction through geologically conditioned petrophysical constraints in joint inversion conditioned petrophysical constraints, *Geophysics*, **82**(6), ID19–ID34.

*   Giraud, J., Lindsay, M., Ogarko, V., Jessell, M., Martin, R. & Pakyuz-Charrier, E., 2019a. Integration of geoscientific uncertainty into geophysical inversion by means of local gradient regularization, *Solid Earth*, **10**(1), 193–210.

*   Giraud, J., Ogarko, V., Lindsay, M., Pakyuz-Charrier, E., Jessell, M. & Martin, R., 2019b. Sensitivity of constrained joint inversions to geological and petrophysical input data uncertainties with posterior geological analysis, *J. geophys. Int.*, **218**(1), 666–688.

*   Giraud, J., Lindsay, M. & Jessell, M., 2021a. Generalization of level-set inversion to an arbitrary number of geologic units in a regularized least-squares framework, *Geophysics*, **86**(4), R623–R637.

*   Giraud, J., Ogarko, V., Martin, R., Jessell, M. & Lindsay, M., 2021b. Structural, petrophysical, and geological constraints in potential field inversion using the tomofast-x v1.0 open-source code, *Geosci. Model Develop.*, **14**(11), 6681–6709.

*   Giraud, J., Caumon, G., Grose, L. & Cupillard, P., 2022. Geometrical inversion coupled with automated geological modelling, in *83rd EAGE Annual Conference & Exhibition*, Vol. **2022**, pp. 1–5, EAGE Publications BV, European Association of Geoscientists & Engineers.

*   Giraud, J., Seillé, H., Lindsay, M.D., Visser, G., Ogarko, V. & Jessell, M.W., 2023. Utilisation of probabilistic magnetotelluric modelling to constrain magnetic data inversion: proof-of-concept and field application, *Solid Earth*, **14**(1), 43–68.

*   Giraud, J., Caumon, G., Grose, L., Ogarko, V. & Cupillard, P., 2024. Integration of automatic implicit geological modelling in deterministic geophysical inversion, *Solid Earth*, **15**(1), 63–89.

*   Gómez-Romeu, J., Masini, E., Tugend, J., Ducoux, M. & Kusznir, N., 2019. Role of rift structural inheritance in orogeny highlighted by the western pyrenees case-study, *Tectonophysics*, **766**, 131–150.

*   Guillen, A., Calcagno, P., Courrioux, G., Joly, A. & Ledru, P., 2008. Geological modelling from field data and geological knowledge, Part II - modelling validation using gravity and magnetic data inversion, *Phys. Earth planet. Inter.*, **171**, 158–169.
Gunnell, Y., Zeyen, H. & Calvet, M., 2008. Geophysical evidence of a missing lithospheric root beneath the eastern pyrenees: Consequences for post-orogenic uplift and associated geomorphic signatures, *Earth planet. Sci. Lett.*, **276**(3), 302–313.

Haber, E. & Oldenburg, D., 1997. Joint inversion: a structural approach, *Inverse Probl.*, **13**(1), 63, doi:10.1088/0266-5611/13/1/006.

Hansen, P., 2001. The L-curve and its use in the numerical treatment of inverse problems, in *Advances in Computational Bioengineering*, Vol. 5, pp. 119–142, ed. Press, W., Computational Inverse Problems in Electrocardiology.

Huyghe, D., Mouthereau, F., Segalen, L. & Furió, M., 2020. Long-term dynamic topographic support during post-orogenic crustal thinning revealed by stable isotope (δ18O) paleo-altimetry in eastern Pyrenees, *Sci. Rep.*, **10**(1), doi:10.1038/s41598-020-58903-w

Jolivet, L., Romagny, A., Gorini, C., Maillard, A., Thinon, I., Couëffé, R., Ducoux, M. & Séranne, M., 2020. Fast dismantling of a mountain belt by mantle flow: Late-orogenic evolution of pyrenees and liguro-provençal rifting, *Tectonophysics*, **776**, 228312, doi:10.1016/j.tecto.2019.228312.

Kadu, A., van Leeuwen, T. & Mulder, W.A., 2017. Salt reconstruction in full-waveform inversion with a parametric level-set method, *IEEE Trans. Comput. Imag.*, **3**(2), 305–315.

Lajaunie, C., Courrioux, G. & Manuel, L., 1997. Foliation fields and 3D cartography in geology: principles of a method based on potential interpolation, *Math. Geol.*, **29**(4), 571–584.

Lefeuvre, N., Truche, L., Donzé, F.-V., Gal, F., Tremosa, J., Fakoury, R.-A., Calassou, S. & Gaucher, E., 2022. Natural hydrogen migration along thrust faults in foothill basins: the north pyrenean frontal thrust case study, *Appl. Geochem.*, **145**, 105396, doi:10.1016/j.apgeochem.2022.105396.

Lelievre, P.G. & Farquharson, C.G., 2016. *Integrated Imaging for Mineral Exploration*, pp. 137–166, c hap. 8, American Geophysical Union (AGU).

Lelievre, P.G., Farquharson, C.G. & Hurich, C.A., 2012. Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration, *Geophysics*, **77**(1), K1–K15.

Lescoutre, R., 2019. *Formation and Reactivation of the Pyrenean-Cantabrian Rift System: Inheritance, Segmentation and Thermal Evolution*, Theses, Université de Strasbourg.

Li, W., Leung, S. & Qian, J., 2014. A level-set adjoint-state method for crosswell transmission-reflection traveltime tomography, *J. geophys. Int.*, **199**(1), 348–367.

Li, W., Lu, W., Qian, J. & Li, Y., 2017. A multiple level-set method for 3d inversion of magnetic data, *Geophysics*, **82**(5), J61–J81.

Li, W., Qian, J. & Li, Y., 2020. Joint inversion of surface and borehole magnetic data: A level-set approach, *Geophysics*, **85**(1), J15–J32.

Li, Y. & Oldenburg, D.W., 2000. Joint inversion of surface and three-component borehole magnetic data, *Geophysics*, **65**(2), 540–552.

Lines, L., 1999. Density and AVO, *Canadian J. Explor. Geophys.*, **35**(1/2), 32–35.

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

Liu, D., Khambampati, A.K. & Du, J., 2018. A parametric level set method for electrical impedance tomography, *IEEE Trans. Med. Imag.*, **37**(2), 451–460.

Ludwig, W., Nafe, J. & Drake, C., 1970. *Seismic Refraction, the Sea*, Vol. 4, ed. Maxwell, A.E., Wiley-Interscience, New-York, USA.

Maag, E. & Li, Y., 2018. Discrete-valued gravity inversion using the guided fuzzy c-means clustering technique, *Geophysics*, **83**(4), G59–G77.

Martin, R., Monteiller, V., Komatitsch, D., Perrouty, S., Jessell, M., Bonvalot, S. & Lindsay, M., 2013. Gravity inversion using wavelet-based compression on parallel hybrid CPU/GPU systems: application to southwest Ghana, *J. geophys. Int.*, **195**(3), 1594–1619.

Martin, R. et al., 2017. A high-order 3-D spectral-element method for the forward modelling and inversion of gravimetric data - Application to the western pyrenees, *J. geophys. Int.*, **209**(1), 406–424.

Martin, R., Ogarko, V., Komatitsch, D. & Jessell, M., 2018. Parallel three-dimensional electrical capacitance data imaging using a nonlinear inversion algorithm and Lp norm-based model regularization, *Measurement*, **128**, 428–445.

Martin, R., Giraud, J., Ogarko, V., Chevrot, S., Beller, S., Gégout, P. & Jessell, M., 2021. Three-dimensional gravity anomaly data inversion in the Pyrenees using compressional seismic velocity model as structural similarity constraints, *J. geophys. Int.*, **225**(2), 1063–1085.

Mistree, F., Lautenschlager, U., Erikstad, S. & Allen, J., 1993. Simulation reduction using the taguchi method. NASA Contractor Report 4542.

Moorkamp, M., Heincke, B., Jegen, M., Hobbs, R.W. & Roberts, A.W., 2016. *Joint Inversion in Hydrocarbon Exploration*, pp. 167–189, chap. 9, American Geophysical Union (AGU).

Mouthereau, F. et al., 2014. Placing limits to shortening evolution in the pyrenees: role of margin architecture and implications for the iberia/europe convergence, *Tectonics*, **33**(12), 2283–2314.

Noble, M.S., Gesret, A. & Belayouni, N., 2014. Accurate 3-D finite difference computation of traveltimes in strongly heterogeneous media, *J. geophys. Int.*, **199**(3), 1572–1585.

Ogarko, V., Giraud, J., Martin, R. & Jessell, M., 2021. Disjoint interval bound constraints using the alternating direction method of multipliers for geologically constrained inversion: application to gravity data, *Geophysics*, **86**(2), G1–G11.

Ogarko, V., Frankcombe, K., Liu, T., Giraud, J., Martin, R. & Jessell, M., 2024. Tomofast-x 2.0: an open-source parallel code for inversion of potential field data with topography using wavelet compression, *Geosci. Model Develop.*, **17**(6), 2325–2345.

Paasche, H. & Tronicke, J., 2007. Cooperative inversion of 2d geophysical data sets: A zonal approach based on fuzzy c-means cluster analysis, *Geophysics*, **72**(3), A35–A39.

Paige, C.C. & Saunders, M.A., 1982. LSQR: an algorithm for sparse linear equations and sparse least squares, *ACM Trans. Math. Softw.*, **8**(1), 43–71.

Pedrera, A., García-Senz, J., Ayala, C., Ruiz-Constán, A., Rodríguez-Fernández, L.R., Robador, A. & González Menéndez, L., 2017. Reconstruction of the exhumed mantle across the north iberian margin by crustal-scale 3-d gravity inversion and geological cross section, *Tectonics*, **36**(12), 3155–3177.

Pedrera, A., García-Senz, J., Ayala, C., Ruiz-Constán, A., Rodríguez-Fernández, L.R., Robador, A. & González Menéndez, L., 2018. Reply to comment by Pedreira et al. on “Reconstruction of the exhumed mantle across the north iberian margin by crustal-scale 3-D gravity inversion and geological cross section”, *Tectonics*, **37**(11), 4346–4356.

Phadke, M., 2021. *Quality Engineering Using Robust Design*, Phadke Associates, Incorporated.

Plazolles, B., Martin, S., Vincent, R. & El Baz, D., 2015. Monte-carlo analysis of object reentry in earth's atmosphere based on taguchi method, in *8th European Symposium on Aerothermodynamics for Space Vehicles*, Lisbonne, Portugal, https://hal.science/hal.01149587.

R.C. Bose, K.B., 1952. Orthogonal array of strength two and three, *Ann. Math. Stat.*, **Vo., 23**, 508–524, doi:10.1214/aoms/1177729331.

Ren, Z. & Kalscheuer, T., 2020. Uncertainty and resolution analysis of 2d and 3d inversion models computed from geophysical electromagnetic data, *Surv. Geophys.*, **41**(1), 47–112.

Roy, K.R., 1990. *A Primer on the Taguchi Method*, Society of Manufacturing Engineers, Michigan.

Sun, J. & Li, Y., 2015. Multidomain petrophysically constrained inversion and geology differentiation using guided fuzzy c-means clustering, *Geophysics*, **80**(4), ID1–ID18.

Sun, J. & Li, Y., 2016. Joint inversion of multiple geophysical and petrophysical data using generalized fuzzy clustering algorithms, *J. geophys. Int.*, **208**(2), 1201–1216.

Taguchi, G., 1987. *System of Experimental Design: Engineering Methods to Optimize Quality and Minimize Costs*, Vol. 1. UNIPUB/Kraus International Publications.

Taguchi, G., Chowdhury, S. & Wu, Y., 2005. *Taguchi's Quality Engineering Handbook*, Wiley.

Torné, M. et al., 2023. Advances in the modeling of the iberian thermal lithosphere and perspectives on deep geothermal studies, *Geotherm. Energy*, **11**(1), 3, doi:10.1186/s40517-023-00246-6.

Wang, Y. et al., 2016. The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic P-waves, *Geology*, **44**(6), 475–478.
Wehr, H., 2017. Modélisation 3D des Pyrénées à partir des données géologiques, gravimétriques et sismiques, PhD thesis, Thèse de doctorat dirigée par Chevrot, Sébastien Sciences de la terre et des planètes solides Toulouse 3 2017.

Wehr, H., Chevrot, S., Courrioux, G. & Guillen, A., 2018. A three-dimensional model of the pyrenees and their foreland basins from geological and gravimetric data, Tectonophysics, 734–735, 16–32.

Zheglova, P., Farquharson, C.G. & Hurich, C.A., 2013. 2-D reconstruction of boundaries with level set inversion of traveltimes, J. geophys. Int., 192(2), 688–698.

Zheglova, P., Lelièvre, P.G. & Farquharson, C.G., 2018. Multiple level-set joint inversion of traveltime and gravity data with application to ore delineation: A synthetic study, Geophysics, 83(1), R13–R30.

## APPENDIX A: INVERSION METHOD CONSTRAINED BY ADMM

Let us describe here our methodology for gravity data inversion under lithological units constraints (Ogarko et al. 2021, 2024; Giraud et al. 2023, 2024). We minimize a cost function $\chi$ that involves the misfit between measured and computed complete Bouguer gravity anomalies, and we regularize the problem by adding prior model and smoothing constraints, as well as geological constraints using density bounded intervals. It is minimized under the constraint that the density anomaly model $m_k$ varies inside a set of possible intervals $B_k$ defined as

$$B_k = \cup_{l=1,L}^{k} [a_{kl}, b_{kl}] \text{ with } 1 \le k \le n \quad (A1)$$

where $B_k$ is the interval binding each of the n model cells. For the kth model cell, $B_k$ is defined as a union of disjointed intervals where $a_{kl}$ and $b_{kl}$ are the lower and upper bounds for kth cell and l is the lithology unit number. The total cost function is thus

$$\chi(m) = \| D_{\text{obs}} - D_{\text{cal}} \|_2^2 + \alpha^2 \| W(\Delta\rho - \Delta\rho_p) \|_2^2 + \alpha_g^2 \| \nabla\Delta\rho \|_2^2, \quad (A2)$$

constrained by the geological bounds $B_k$ such that each $m_k = \Delta\rho_k$ model cell verifies $m_k \in B_k$. $\Delta\rho$ is the density anomaly to be retrieved by inversion and $\Delta\rho_p$ the prior density anomaly relative to a reference density model $\rho_{\text{ref}}$ defined as follows:

$$\begin{align*}
\rho_{\text{ref}} &= \rho_{\text{crust}} = 2670 \, \text{kg m}^{-3}, & \text{for a depth } \le 29 \, \text{km} \\
\rho_{\text{ref}} &= \rho_{\text{upper mantle}} = 3260 \, \text{kg m}^{3}, & \text{for a depth } > 29 \, \text{km} \tag{A3}
\end{align*}$$

$D_{\text{obs}}$ and $D_{\text{cal}}$ are the measured and theoretical Bouguer gravity data, and $\alpha$ and $\alpha_g$ the regularization weights applied to the density model and the model gradients, respectively, to simultaneously smooth models and take into account the gradients at interfaces.

The depth weighting operator $W$ is applied on the model by computing for each cell volume $\Delta V_k$, and it is defined by the distance between a cell and all observation locations. For the density anomaly in a cell $k$ of the computational mesh, this matrix operator $W$ is defined as in Li & Oldenburg (2000):

$$W_{kk} = \left\{ \sum_{i=1}^{N_{\text{obs}}} \left[ \int_{\Delta V_k} \frac{dv}{(R_{ik} + R_0)^2} \right]^2 \right\}^{1/4}, \quad k = 1, \dots, n, \quad (A4)$$

where $R_{ik}$ is the distance between the i'th observation and any point within $\Delta V_k$. The parameter $R_0$ is a small constant, chosen here to be 0.1m.

Following Boyd et al. (2011) and the related formulations more adapted to our geophysical application as in Ogarko et al. (2021),

this problem is solved iteratively with the ADMM method by alternating the updates of the model $m$ and temporary variables $z$ and $u$ as:

$$\begin{align}
m^i &= \operatorname{argmin}(\chi(m) + \alpha_{\text{admm}}^2 \| m - z^{i-1} + u^{i-1} \|_{L^2}^2) \nonumber \\
z^i &= P_B(m^i + u^{i-1}) \tag{A5} \\
u^i &= u^{i-1} + m^i - z^i \nonumber
\end{align}$$

where $\alpha_{\text{admm}}$ is a penalty parameter. The bounds constraint $z$ is obtained as the projection of a perturbation $p = m + u$ of the model $m$ on the different interval sets $B$ such that

$$P_B(p) = [P_{B1}(p_1), \dots, P_{Bk}(p_k), \dots, P_{Bn}(p_n)] \quad (A6)$$

with

$$P_{Bk}(p_k) = \operatorname{argmin} \| p_k - x \|_2 \text{ for } x \in B_k \quad (A7)$$

In our Tomofast-x parallel inversion code, we minimize equation (A2) using the LSQR algorithm of Paige & Saunders (1982), which is implemented using a parallelized sparse matrix solver. At each inversion cycle (i.e. outer loop of LSQR inversion), the model is updated using an approximate solution using a fixed number of iterations in the inner loop of the LSQR algorithm (Martin et al. 2013, 2018, 2021; Giraud et al. 2019a, b; Ogarko et al. 2021, 2024). A least-squares algorithm is used here to solve the matrix system as follows. Let us denote $m$ the model vector containing the density anomalies $\Delta\rho$. This model vector is obtained by solving the following system of equations:

$$\begin{bmatrix} SW^{-1} \\ \alpha I_d \\ \alpha_g W^{-1} \nabla \end{bmatrix} [Wm] = \begin{bmatrix} \Delta g_{\text{obs}} \\ \alpha W m_p \\ 0 \end{bmatrix}, \quad \text{under the constraint that } m_k \in B_k \quad (A8)$$

where $m_p = \Delta\rho_p = \rho_p - \rho_{\text{ref}}$ is the a priori density anomaly relative to a reference model $\rho_{\text{ref}}$ defined as in eq. (A3) and the gravity kernel $S$ is depth-weighted by the operator $W$ given in eq. (A4). The solution of eq. (A8) is obtained by an iterative LSQR algorithm. At each i'th LSQR inversion cycle we solve the following linearized version of the system of eqs (A8)

$$\begin{align*}
& \text{For } 1 \le i \le \text{Niter}_{\text{max}} \\
& \quad \begin{bmatrix} SW^{-1} \\ \alpha I_d \\ \alpha_g W^{-1} \nabla \\ \alpha_{\text{admm}} I_d \end{bmatrix} [\Delta \bar{m}^i] \\
& \quad = \begin{bmatrix} \Delta g_{\text{obs}} - \Delta g(m^{i-1}) \\ -\alpha W(m^{i-1} - m_p) \\ -\alpha_g \nabla m^{i-1} \\ -\alpha_{\text{admm}} W(m^{i-1} - z^{i-1} + u^{i-1}) \end{bmatrix} \\
& \quad \Delta m^i = W^{-1} \Delta \bar{m}^i \\
& \quad m^i = m^{i-1} + \Delta m^i \\
& \quad \text{End for} \\
& \quad , \tag{A9}
\end{align*}$$

where $\text{Niter}_{\text{max}}$ is the maximum number of outer loop inversion cycles, $u^0 = 0$, $z^0 = 0$, and the starting model $m^0 = \Delta\rho^0 = 0$. The updated ADMM variables $z^i$ and $u^i$ are calculated using eqs (A6)–(A9). Spatial gradients are computed using forward finite-difference integration scheme. In order to make the whole algorithm converge under the ADMM constraint it is crucial to take a small initial $\alpha_{\text{admm}}^0$ value and make it increase during the inversion process, this way both data misfit and ADMM cost functions are converging with similar/equilibrate values.
<figure><img src="image_26.png" /></figure>

**Figure A1.** N-S sections of the difference between the inverted and prior models for a 3150 kg m⁻³ IUM prior density (2 km resolution along the first four south-north seismic PYROPE vertical profiles). Thick dotted line represents the crust/mantle boundary. Units of both abscissa and vertical coordinates are in kilometres. Negative density values show that the prior model overestimates density with respect to the inverted model. Positive density values show that the prior model underestimates density with respect to the inverted model. We applied a conservative padding value of 20 km (in N-S profiles) at the beginning and the end of each profile (white shaded area) to avoid edge effects. Note that the E-E' profile has been discarded because of too many artefacts due to poor data gravity data and possible effects of the Gulf of Lion. A-H: Aston-Hospitalet; AZ-3S: Arize-Trois Seigneurs; LA: Labourd anomaly; NPF: North Pyrenean Fault; NFPT: North Pyrenean Frontal Thrust; NPZ: North Pyrenean Zone; SG: Saint Gaudens anomaly and SPFT: South Pyrenean Frontal Thrust.
<figure><img src="image_27.png" /><figcaption>**Figure A2.** E-W sections of the difference between the inverted model and the prior model for a 3150 kg m⁻³ IUM prior density. Thick dotted line represents the crust/mantle boundary. Negative density values show that the prior model overestimates density with respect to the inverted model. Positive density values show that the prior model underestimates density with respect to the inverted model. We applied a conservative padding value of 40 km (in E-W profiles) at the beginning and the end of each profile (white shaded area) to avoid edge effects. Note that the E-E' profile has been discarded because of too many artefacts due to poor data gravity data and possible effects of the Gulf of Lion. A-H: Aston-Hospitalet; AZ-3S: Arize-Trois Seigneurs; LA: Labourd anomaly; NPF: North Pyrenean Fault; NFPT: North Pyrenean Frontal Thrust; NPZ: North Pyrenean Zone; SG: Saint Gaudens anomaly-and SPFT: South Pyrenean Frontal Thrust.</figcaption></figure>