Basement and cover architecture in the Central Pyrenees constrained by gravity data

A new gravity survey (1164 gravity stations and 180 samples for density analysis) combined with two new geological cross sections has been carried out in a sector of the Central Pyrenees in order to improve the characterization of basement and cover architecture. From North to South, the study area comprises the southern half of the Axial Zone and the northernmost part of the South-Pyrenean Zone. New gravity data were combined with previous existing databases to obtain the Bouguer and residual anomaly maps of the study area. The two cross sections, oriented NNE–SSW, were built from field data and previous surficial and subsurface data and cross the La Maladeta plutonic complex. The residual anomaly map shows values ranging from −18 to 16 mGal and anomalies mainly oriented N120E. The two 2.5D modelled cross sections show similar observed gravity curves coinciding with similar interpreted structural architecture. Data show a gravity high oriented N120E coinciding with the Orri basement thrust sheet and an important gravity depression, with the same orientation, coinciding with the leading edge at depth of the Rialp basement thrust sheet and interpreted as linked to a large subsurface accumulation of Triassic evaporites. The volume at depth of the La Maladeta and Arties granites has been constrained through gravity modelling. This work highlights that the combination of structural geology and gravity modelling can help to determine the structural architecture of an orogen and localize accumulations of evaporites at depth.


Introduction
The Pyrenees (Northeast of Spain and South of France) have attracted multitude of research during the past decades due to their excellent outcrops generating an extensive knowledge of its surface geology and map-based reconstructions. Regarding its subsurface geology, the ECORS-Pyrenees deep seismic reflection profile acquired during the 80s provided the first image of the structure and reflectivity of the crust across the Central Pyrenees. The results show an asymmetric doubly verging thrust wedge that formed above the subduction of the Iberian crust under the European plate (Chokroune and ECORS Team 1989;Roure et al. 1989;Torné et al. 1989;Muñoz 1992). Since then, numerous geophysical studies have been performed, with different interpretations leading to quite contrasting crustal models (e.g. Chevrot et al. 2014Chevrot et al. , 2018Wehr et al. 2018;Pedrera et al. 2017;Teixell et al. 2018;García-Senz et al. 2019; and references therein).
In the central part of the chain, the internal architecture of the orogen consists of an antiformal stack of basementinvolved units dominated by Paleozoic rocks (Axial Zone) and two fold-and-thrust belts involving Mesozoic and Cenozoic units flanking the Axial Zone to the north and south (e.g., Muñoz 1992; Fig. 1). The Axial Zone is characterized by the presence of several Late Variscan igneous bodies emplaced in the Paleozoic metasedimentary succession (Zwart 1986; Barnolas et al. 1996). In large areas to the south of the Axial Zone, the Mesozoic and Cenozoic cover was decoupled from the basement along the Triassic evaporites, which act as the regional décollement. The subsurface structure of the Mesozoic and Cenozoic cover has been mainly constrained from the interpretation of seismic data and wells acquired during the 60s, 70s and 80s (e.g., Muñoz 1992;Vergés 1993;Teixell and Muñoz 2000;Muñoz et al. 2018). However, in the Axial Zone and leading edges, subsurface data are very scarce and the geometry of structures at depth is more speculative.
An efficient and cost-effective way to better constrain the subsurface structure in areas without seismic reflection and/ or wells is the use of potential field geophysics. Gravimetry is a geophysical technique sensitive to lateral density variations thus providing a constraint to reconstruct geometries at depth (e.g. Goleby et al. 1989;Ayala et al. 2019). Since Torné et al. (1989), several works including potential field modelling have been published trying to image the deep structure of the Pyrenees down to the lithospheric mantle (e.g., Adam 1993; Casas et al. 1997;Ledo et al. 2000;Vacher and Souriau 2001;Jammes et al. 2010;Pedrera et al. 2017Pedrera et al. , 2018Chevrot et al. 2018;Wehr et al. 2018;García-Senz et al. 2019) or focusing on the upper crustal structures covering more specific zones in the southern Pyrenees (e.g. Santolaria et al. 2016;Calvín et al. 2018).
Considering the density contrast between host rocks and granites in the basement and between evaporites and other sedimentary rocks in the cover along the Central Pyrenees, in this work we use a workflow that combines the construction of two geological cross sections and their corresponding gravimetric models (see also Santolaria et al. 2016Santolaria et al. , 2020Izquierdo-Llavall et al. 2019). The aim of this work was to investigate the basement and cover architecture in the  Muñoz (1992 Central Pyrenees and to constrain the geometry at depth and the volumetric distribution of the Triassic evaporites and the La Maladeta batholith. Furthermore, we intend to test if the overall proposed geometry of the Central Pyrenees is consistent with the gravity data. The target area is a sector of the Central Pyrenees encompassing the southern half of the Axial Zone and the northernmost part of the South-Pyrenean Zone (Fig. 1). This area is key for understanding the role played by the Triassic evaporites on the mechanical decoupling between the basement and cover rocks.

Geological setting
The Pyrenees constitute an Alpine range that resulted from the collision between the Iberian and European plates from the Late Cretaceous to the Miocene (e.g. Muñoz 1992). It is an asymmetric doubly vergent orogen with a mainly southward vergence (e.g. Muñoz 1992) (Fig. 1b). Classically the range has been divided into the following three WNW-ESE trending structural zones (Mattauer 1968;Fig. 1a): (i) The North-Pyrenean Zone, (ii) the South-Pyrenean Zone, and (iii) the Axial Zone at the core of the range. The North-Pyrenean Zone mainly consists of Mesozoic sediments overlain by syncompressional Paleogene sediments (Déramond et al. 1988) characterized both by a mechanical coupling between Paleozoic basement and cover rocks (thick-skinned tectonics) in some sectors and a total decoupling resulting from the presence of a thick detachment level (clays and gypsum) at the Upper Triassic in other areas. The South-Pyrenean Zone is characterized by a thicker Paleogene syntectonic succession (e.g. Séguret 1972;Muñoz 1992), a clear southward vergence and a more generalized decoupled deformation between the cover and the basement. Deformation of the Pyrenean double-wedge progressed outward in a piggyback fashion, although synchronous inward internal deformation has been also documented (e.g. Muñoz et al. 1986;Vergés and Muñoz 1990). The Southern Central Pyrenees are characterized by a central antiformal stack of several southward-facing basement-involved thrust sheets (i.e. the Axial Zone) and an imbricate stack involving cover rocks (i.e. the South-Pyrenean Zone) (e.g. Muñoz 1992;Beaumont et al. 2000). This structural architecture responds to the interaction between a mid-crustal detachment effective for the basement-involved thrust sheets and the Triassic evaporites that decoupled deformation between the Mesozoic succession and the basement in the South-Pyrenean Zone (e.g. Muñoz 1992;Beaumont et al. 2000). The Axial Zone is composed of several basement sheets whose number, name and distribution vary along-strike and from an author to another (e.g. Parish 1984;Williams 1985;Poblet 1991;Muñoz 1992;García-Sansegundo 1992;Teixell 1996;Martínez-Peña and Casas-Sainz 2003;Labaume et al 2016;Izquierdo-Llavall et al. 2018;García-Senz et al. 2019;Espurt et al. 2019). In this work, in the study area, we have considered the following main basement thrust sheets differentiated by Muñoz (1992): the Nogueres, Orri and Rialp thrust sheets (Fig. 1b) and the Ribagorza thrust sheet (Muñoz et al. 2018). They were progressively deformed by underthrusting of the lower and younger units (e.g. Muñoz 1992). The uppermost thrust sheet (the Nogueres Unit) was forward rotated and tilted in the southern limb of the antiformal stack showing downward facing structures (Séguret 1972;Muñoz 1992). Rocks belonging to basement thrust sheets were already deformed during the Variscan Orogeny (Middle-Late Carboniferous), whose main structural features were partly reactivated during the Alpine compression (e.g. Poblet 1991;García-Sansegundo 1996;Gil-Peña 2004). During the Alpine shortening, the Triassic evaporites acted as the upper décollement for the Axial Zone antiformal stack (i.e. roof thrust) (e.g. Muñoz 1992;Saura and Teixell 2006). Internally, the Silurian black shales also acted as a décollement both in the Variscan and Alpine Orogenies (Matte 1969;García-Sansegundo 1990;García-Sansegundo et al. 2011).
South of the Axial Zone in the Central Pyrenees, the Mesozoic-Cenozoic cover forms a large imbricate thrust system known as South Pyrenean Central Unit (SPCU; Séguret 1972) (Fig. 1a). The SPCU consists of three main thrust units (from top to base and in order of emplacement, Bóixols, Montsec and Sierras Marginales thrust sheets) (Fig. 1b) emplaced in a piggyback thrust sequence from the Late Cretaceous to the Oligocene. The study area only encompasses the northernmost unit, the Bóixols thrust sheet (Fig. 1a, b). The Bóixols unit resulted from the inversion of one of the main Lower Cretaceous extensional basins of the Pyrenees, the Organyà basin, from Late Santonian to Maastrichtian (Garrido-Megías 1973;Simó 1986;Bond and McClay 1995;García-Senz 2002). Between the Bóixols thrust sheet and the Axial Zone, it appears a sector, named Ribagorza basin (Fig. 2), characterized by several large Triassic salt diapirs (e.g. Aulet diapir) and a complex fault network (Teixell and Muñoz 2000;García-Senz 2002;Saura et al. 2016). The Ribagorza basin consists of a Cretaceous basin inverted during the Pyrenean orogeny (Saura et al. 2016). The Triassic evaporites, at the base of the Mesozoic succession, played a major role on the structural architecture both during the Mesozoic rifting and its posterior tectonic inversion (e.g. García-Senz 2002;Mencos et al. 2015;Muñoz et al. 2018). Southwards of the study area, Triassic outcrops are scarce and this unit re-appears in the large exposures located within the Sierras Marginales thrust sheet, the youngest unit of the SPCU.

Undiferenciated fault
Unconformities N o g u e ra R ib a g o rz a n a

Stratigraphy
The study area comprises a large variety of rocks both in age and lithology. In general terms, the oldest rocks (i.e. basement) crop out in the Axial Zone, located in the northern half of the study area, and the cover and younger sequences appear in its southern half (Fig. 2). The stratigraphy of the study area is described considering the division on geological units done for the gravity modelling based on their density values (see sixth section and Fig. 2). The oldest rocks, Cambro-Ordovician in age, are a monotonous alternation of quartzites and slates with some intercalations of limestone and microconglomerate. Above, the Upper Ordovician sequence is mainly composed of siliciclastic rocks that lie unconformably upon the Cambro-Ordovician rocks (García-Sansegundo et al. 2004. In this work the Upper Ordovician sequence has been included within the Cambro-Ordovician unit. The Silurian rocks are represented by black shales that often acted as a décollement, conditioning deformation (Matte 1969;Poblet 1991;García-Sansegundo 1990, 1996. The Devonian units are mainly formed by slates and limestones (Zwart 1979). The youngest rocks of the variscan succession consist mainly of a siliciclastic package (Culm facies; Devolvé 1987). Upper Pennsylvanian and Permian rocks lie unconformably upon deeply eroded older rocks. The Upper Pennsylvanian units are mainly formed by andesitic lava flows and pyroclastic units with some intercalated coal beds. The Permian succession is made of sandstone, shales and microconglomerates (e.g. Mey et al. 1968;Gisbert 1983). In the study area several igneous bodies crop out (the La Maladeta plutonic complex, the Arties stock, the Marimanha pluton, the Bossòst leucogranite and the Lys-Caillaouas pluton) (Fig. 2). They intruded variscan rocks and some of them, including the La Maladeta body have been dated to be emplaced during the Late Carboniferous-Early Permian (Evans 1993;Solé et al. 1997;Esteban et al. 2015;Mezger and Gerdes 2016;López-Sánchez et al. 2019). The Triassic succession shows the following typical Germanic facies: Buntsandstein red beds, Muschelkalk dolostones and limestones and Keuper evaporites and shales. Hypabyssal rocks (dolerites) are relatively common within the Upper Triassic rocks (Lago et al. 2000). In this work, Upper Permian rocks and the Buntsandstein red beds have been grouped together and the Muschelkalk and Keuper are separated in a different unit.
The Jurassic and Lower Cretaceous sequences were deposited during a rifting period related to the opening of the Bay of Biscay and North Atlantic Ocean (e.g. Roca et al. 2011) and are represented by marine limestones and marls (Aurell and Meléndez 2002;García-Senz 2002). The Upper Cretaceous sedimentary deposits consist of a succession of limestones, calcarenites and marls (Simó 2004) that registered an evolution from marine to transitional and finally continental environments (Puigdefàbregas and Souquet 1986;Puigdefàbregas et al. 1992). The limit between the Mesozoic and Cenozoic successions is located within a continental unit (Garumnian facies) consisting of red clays, sandstones and lacustrine limestones (Pujalte and Schmitz 2005). In this work we have merged this unit with the Lower Paleogene. This unit also comprises limestones, sandstones and marls belonging to the Graus-Tremp basin . The youngest deposits considered in the study area, here included in the Upper Paleogene unit, are unconformable continental conglomerates deposited during Late Lutetian to Late Oligocene times (Beamud et al. 2011).

Methodology
Two parallel geological cross sections oriented NNE-SSW and separated approximately 12 km were built to infer the structural architecture of the study area (MAL and AG1; Fig. 2). They run perpendicular to the trend of the main geological structures. They were constructed using field data, surface geological data available from geological maps (Rosell 1994 , as well as data from four nearby oil exploration wells (Tamurcia-1 -Ta-1-, Cajigar-1 -Cj-1-, Comiols-1 and Isona-1, Lanaja 1987). The two first boreholes ( Fig. 2) reached the décollement level located on the Triassic evaporites and the other two (located to the South of the study area, out of the frame of Fig. 2) reach the autochthonous Eocene rocks allowing to know the minimum depth to the top of the basement (Teixell and Muñoz 2000). The geometry at depth interpreted by several geological cross sections already published by different authors was considered (references in Fig. 2).
Gravity data comes from the Institut Cartogràfic de Catalunya (CAT AGR AV database 2010) and several gravity surveys carried out in the frame of GeoPiri3D project in 2018 and 2019  to improve the spatial resolution of the data. After merging the datasets (CAT AGR AV, Ayala et al. 2020 and the newly acquired data), approximately 1 station every km 2 is available surrounding La Maladeta and Arties granites, except in some areas above 2000 m a. s. l. that were inaccessible by car (Fig. 3a). To have a homogeneous gravity data (Fig. 3b), we have calculated the complete Bouguer anomaly of the new data with the same parameters of the following databases: GRS80 geodetic system formulae, orthometric heights and a density reduction of 2.67 g/cm 3 . The following corrections were applied to the theoretical gravity to obtain the Bouguer anomaly of the new data: free-air correction, Bouguer slab correction and terrain correction. The terrain correction was applied up to 167 km using the Oasis Montaj ® terrain correction module that is based on Kane (1962) and Nagy (1966) algorithms, using a combined 5 m DTM from Spain, France and Andorra that were regridded to 30 × 30 m. as local grid and a combined grid from SRTM database and bathymetry data from GEBCO database regridded to 100 × 100 m as regional grid. To carry out the analysis and interpretation of the Bouguer anomaly, we have gridded the data with a grid spacing of 1000 m using the minimum curvature method.
The cross sections were improved through the 2D-2.5D gravity modelling in a feed-back process between the geological and gravimetric modelling until the gravimetric response of the models fitted the observed gravity data. Based on the geological surface data and taking into account the aim of this study, we decided to limit the lateral extent only of the outcropping plutonic rocks that have clear boundaries with the host rocks, so only the La Maladeta granite has been modelled in 2.5D while the host rocks and all other units have been modelled in 2D. Along MAL cross-section, La Maladeta granite has a lateral extension of 8.57 km to the W and 32.64 km to the E. Along AG1 cross-section, the extension of La Maladeta granite is 21.18 km to the W and 19.57 km to the E. The modelling was performed using the GM-SYS module of the Oasis Montaj ® software by Seequent based on the Talwani et al. (1959) and Won and Bevis (1987) algorithms. The models were extended far enough to avoid edge effects. Since the target was the upper crust, we have used the residual Bouguer anomaly as observable.

Bouguer and residual anomaly maps
The new Bouguer anomaly map used in this paper (Fig. 3b) is characterized by a long wavelength elongated minimum of roughly E-W direction that occupies the central part of the study area with an amplitude of ca. −100 mGal. This minimum, associated with the crustal root of the Pyrenees (e. g. Torné et al. 2015), seems to continue towards the W and closes towards the E (outside of the study area), presenting on the edges a positive gradient reaching Bouguer anomaly values up to −65 mGal. Different relative maxima and minima of medium and short wavelengths with variable amplitudes are superimposed on this minimum and can be associated with shallower structures.
La Maladeta granite is a massive batholith that crops out in the central Pyrenees. The Bouguer anomaly over the La Maladeta granite outcrop shows variations in amplitude of only 12 mGal (Fig. 3b). There are other smaller granites in this area whose outcrops are associated with a gravimetric minimum. One of them is the Arties granite (Fig. 3b), to the north of the La Maladeta granite, that is characterized by a relative minimum of −100 to −108 mGal within the northern gradient zone. The relative minimum located to the SW of La Maladeta has been interpreted as the gravimetric response of Triassic evaporitic accumulations far greater than the mapped outcrops.
The residual anomaly (Fig. 3a), that reflects the structures of the upper crust which are the targets of our study, was calculated by subtracting from the Bouguer anomaly, a regional anomaly that we assumed to correspond to a third degree polynomial (Fig. 3c). Focusing on the study area, La Maladeta granite shows a gravimetric zonation with small variations in its amplitude from one zone to the next (from 6 to −10 mGal) (Fig. 3a), consistent with lateral changes in its composition, predominantly granodioritic (Debon et al. 1996). The anomaly over the Arties granite has a minimum that reaches ca. −9 mGal, located in its westernmost part. From there, the anomaly increases steadily towards the borders with a rather gentle gradient reaching values up to −2 mGal. The relative minimum with WNW-ESE orientation associated with the Triassic evaporites is one of the most prominent features of the residual map with a total amplitude of ca. 15 mGal (Fig. 3a).

Density measurements
The characterization of the rock density is fundamental in gravity surveys to constrain the modelling thus reducing the uncertainties in the modelled geometry (e. g. McCulloh 1965;Santolaria et al. 2016;Izquierdo-Llavall et al. 2019). The densities were measured from rock samples of the different lithologies represented in the two cross sections. A total of 180 samples were collected from non-weathered outcrops which represent the whole range of geological units cropping out within the study Fig. 3 a Residual Bouguer anomaly map and distribution of the gravity stations by source: red circles indicate data from the CAT AGR AV database (ICGC 2010); black circles correspond to data from Ayala et al. (2020) complemented in this study. Black lines indicate the location of the modelled cross sections. The contour lines correspond to La Maladeta (1), Arties (2), Marimanha (3), Lys-Caillaouas (4) and Bossost (5) granites. b Bouguer anomaly map of the Central Pyrenees modified from Ayala et al. 2020 showing location of the study area (white rectangle). c Regional anomaly of the study area ◂ area (many of them compiled in Ayala et al. 2020). Bulk density was measured using the double weighting method (standard UNE-EN 1936:2007, by weighing irregular samples weighting between 1000 and 2000 g in air and water (Archimedes principle). Taking into account these density values and previously published density data in the Pyrenees (Santolaria et al. 2016(Santolaria et al. , 2020, for the Cenozoic units), we divided the stratigraphic succession of the modelled cross sections into 13 geological units (see Table 1). Individual sample density values range from 1.93 g/cm 3 (Middle-Upper Triassic evaporites) to 3.03 g/ cm 3 (Devonian metamorphosed slates and limestones).
The density values used in the gravity modelling are plotted in the histograms of Fig. 4 where the range, mean and standard deviation of each geological unit have also been displayed together with a graph showing the percentage of the different lithologies (Fig. 4, Table 1). Available subsurface density data of Jurassic (2.7 g/cm 3 ) and Cretaceous (2.65-2.7 g/cm 3 ) limestone from borehole Cajigar 1 (Lanaja 1987;Santolaria et al. 2016) (Cj-1 in Fig. 2) confirm that the subsurface density values are similar to the outcropping samples. Additionally, the petrophysical data used in the modelling are similar, with very few exceptions, to the data obtained in similar studies carried out in other areas of the Northeast of Iberian Peninsula such as the Iberian Ranges (ca. 1500 points) (Pueyo et al. 2016).
To build the 2D-2.5D gravity models and to fit the gravity anomalies over the models, we choose density values as close as possible to the mean values of each unit and always within the range up to ± 0.23 g/cm 3 std dev (see Table 1).

Geological structure based on previous studies
The overall structure observed from the constructed cross sections in the study area shows that almost all structures have a southern vergence. The geological structure and internal architecture of both cross sections are similar in terms of the position, number and characteristics of major structures and thrust sheets. Both cross sections display the same inhomogeneous coverage by geophysical data. In the South, the geometry at depth of the Bóixols thrust sheet has been well established by several authors based on seismic data and wells (Teixell and Muñoz 2000;García-Senz 2002;Mencos et al. 2015;Saura et al. 2016;Muñoz et al. 2018). However, in the Axial Zone only seismic data obtained from the ECORS-Pyrenees deep seismic reflection profile are available and most structural interpretations are mainly based on field data (e.g. Poblet 1991;García-Sansegundo 1992).
With respect to the Axial Zone, in the study area, its structure is defined by north-dipping thrusts with intermediate to steep dips that cut across, or re-activate Variscan structures (Gil-Peña 2004;Cochelin et al. 2017) (Fig. 5). Alpine structures involving the Paleozoic rocks are in general parallel to the Variscan fabric showing E-W to ESE-WNW trend and northward dip (Poblet 1991;García-Sansegundo 1992;Gil-Peña 2004;Gutierrez-Medina 2007). In the study area, the Axial Zone is formed by the following four basement thrust sheets that from top to bottom in the pile: Nogueres, Orri, Rialp and Ribagorza thrust sheets (Muñoz et al. 2018) (Fig. 5). The northern basement thrust sheets show steeper dips than the southern ones (Muñoz 1992;Teixell et al. 2018).
The upper basement thrust sheet in the study area is the Nogueres thrust sheet (Figs. 1, 5). The Nogueres Unit was emplaced during the Late Cretaceous-Middle Eocene and was progressively deformed by the stacking of the lower and younger units (Muñoz 1992;Vergés et al. 1995;Saura 2004). As a consequence of this deformation, the structures of this thrust sheet were tilted in the southern limb of the antiformal stack showing downward facing folds (the socalled têtes plongeantes; Séguret 1972) that characterize the Nogueres Zone (Muñoz 1992) (Fig. 5). In the study area, the Nogueres Zone displays two small thrust slices involving Paleozoic (Devonian rocks, Upper Pennsylvanian volcanic rocks and Permian red beds) and Triassic rocks (red beds and dolostones) (Muñoz 1992;Saura 2004  structures located to the North. Muñoz (1992) proposed that the root of the allochtonous Nogueres Zone corresponds to the Gavarnie thrust (Figs. 2, 5), implying a large displacement along this thrust. This interpretation is discussed by Soler et al. (1998), Laumonier (2015 and Cochelin et al. (2017), who interpret that the Gavarnie thrust had a much lower displacement. Recently, Teixell et al. (2018) proposed that the root the Nogueres Zone is located in the North Pyrenean Fault, interpreted as a south directed thrust steepened by later deformation. Nevertheless, an important drawback for this hypothesis is that at present the hanging wall of this fault is only constituted by post Paleozoic rocks. Thus, the root of the Nogueres thrust sheet is nowadays still a matter of debate. The Gavarnie thrust crops out in the northern part of both cross sections (Fig. 5). It separates Devonian-Silurian slates and limestones of the Aran Valley syncline in the hanging wall from Cambro-Ordovician slates and quartzites in the footwall (Fig. 5), representing a stratigraphic omission (Soler et al. 1998;García-Sansegundo and Ramírez Merino 2013). In cross section MAL (Fig. 5) highly deformed Devonian rocks affected by contact metamorphism (García-Sansegundo 1992) crop out in the hanging wall of the Gavarnie thrust, whereas in cross section AG1 the Arties granite appears (Fig. 5). The Orri thrust sheet (Muñoz 1992) is well exposed in the study area containing large outcrops of Cambro-Ordovician to Lower Carbonifeous rocks unconformably overlain by discontinuous outcrops of Lower-Middle Triassic red beds and intruded by the La Maladeta granite. The Orri thrust sheet groups three main south-directed thrusts, the Orri s.s. (outcropping to the east of the study area), Bono and Llavorsí-Senet thrusts (Poblet 1991;Muñoz 1992;García-Sansegundo 1996;Saura and Teixell 2006). The Bono thrust represents an Alpine structure and crops out in cross section MAL (Fig. 5). It superimposes Lower Devonian rocks over Lower-Middle Devonian rocks and Lower Triassic red beds (Figs. 2, 5) (Gutierrez-Medina 2007;Izquierdo-Llaval et al. 2013). Eastward, in the cross section AG1, this thrust is not recognized at surface (García-Senz and Ramirez-Merino 2009). The Llavorsí-Senet thrust crops out in cross section AG1 to the South of the La Maladeta granite (Figs. 2, 5). It is variscan in age and reactivated during late-variscan and alpine phases as demonstrated by the cross-cutting relations with the contact metamorphic aureole of the La Maladeta granite (Poblet 1991). The northern sector of the Orri thrust sheet, between the Gavarnie thrust and the La Maladeta granite, involves Lower Carboniferous rocks characterized mainly by a siliciclastic sequence "Culm facies" affected by variscan deformation (García-Sansegundo 1992) (Figs. 2,  5). The La Maladeta granite, whose gravimetric characterization will be detailed in the next section, is much larger in cross section AG1 (Fig. 5). To the South of the La Maladeta granite and still located in the Orri thrust sheet, the Baliera-Flamisell Unit (Mey 1967;Mey et al. 1968) can be differentiated. This unit is mainly composed of outcropping Silurian and Devonian rocks deformed by two Variscan deformation phases and some Alpine thrusts (e.g. the Bono thrust) (Poblet 1991;Gil-Peña and Barnolas 2001). The most conspicuous structures in this unit were formed in the second phase of Variscan deformation and are characterized by east-west-trending upright to inclined folds and thrust faults. These thrust emanate from detachment levels located in different stratigraphic positions, both within the Silurian shales and within the Cambro-Ordovician rocks (García-Sansegundo 1990, 1996Poblet 1991;García-Sansegundo et al. 2011).
The younger basement thrust sheets, the Rialp and Ribagorza thrust sheets, do not crop out in the study area. The leading edges at depth of these younger basement thrust sheets are located below the cover thrust sheets (Teixell and Muñoz 2000;Muñoz et al. 2018) (Fig. 5). These positions have been constrained taking into account (i) the depth of the autochthonous Paleogene succession drilled by the Comiols-1 and Isona-1 wells (located to the South of the study area), (ii) the thickness variations of the Mesozoic units and (iii) the identification of seismic reflectors related to the Bóixols unit (Teixell and Muñoz 2000;Muñoz et al. 2018). In this work the leading edges of both basement thrust sheets have been refined based on gravity data.
The Mesozoic succession of the study area is characterized by a widespread layer of Triassic evaporites at the base (Muñoz 1992;Teixell and Muñoz 2000;Muñoz et al. 2018). Cross sections go across the inverted Ribagorza basin (Saura et al. 2016) and the Bóixols thrust sheet, where a complex passive backthrust develops near its boundary with the Axial Zone (Morreres backthrust; Muñoz 1992). Triassic diapirs in the Ribagorza basin show large outcrops and a complex inner structure in which basic hypabyssal rocks (dolerites; locally called "ofites") are relatively common (Lago San José and Pocoví-Juan 1980;Saura et al. 2016), but despite their high density do not show a gravimetric signature due to their small volume. The thickness of the Mesozoic succession increases notably towards the South in these cover thrust sheets (Teixell and Muñoz 2000;García-Senz 2002) (Fig. 5). The Bóixols thrust sheet is bounded by the southdirected Bóixols thrust, which corresponds to a blind structure in the study area. The southern termination of both cross sections displays Paleogene rocks lying unconformably above strata of the previous units (Fig. 5).

2.5D gravity modelling
The geological cross-sections MAL and AG1 were further constrained by means of 2.5D gravity modelling. The RMS (Root Mean Square) of the final models is 1.1 mGal for the MAL profile and 1.6 mGal for the AG1 profile.
The observed gravity data (residual gravity) display a similar pattern along both cross sections (upper part of Fig. 6), which matches well with our interpretation of a similar structural architecture for both cross sections. The observed gravity curves show a relative gravity maxima (up to 2.4 mGal) coinciding with the Orri basement thrust sheet and an important relative gravity minima (−18 mGal in cross section MAL and −17 mGal in cross section AG1) coinciding with the southern termination at depth of the Rialp basement thrust (Fig. 5).
The wide gravity high of the Orri basement thrust sheet shows a roughly flat geometry located between the southern end of the La Maladeta granite and the beginning of the Nogueres Zone that coincide with the Baliera-Flamisell Unit (Figs. 5, 6). This zone of the profile belongs to a N120E elongated anomaly that ranges from −2 to 3.8 mGal, as can be seen on the residual anomaly map shown in Fig. 3a. It correlates with large outcrops of Devonian rocks deformed by numerous south-vergent folds and thrusts and at depth with thinner Devonian and Silurian layers and shallow Cambro-Ordovician sequences (ca. 1000 m below the surface), as deduced from the gravity model. Towards the North, the gravity high shows a gradual decrease in the observed gravity data from the southern end of the La Maladeta granite to the North (Figs. 2, 5, 6). The negative (absolute) values are lower in cross section MAL with respect to cross section AG1 (Fig. 6). This negative decrease is compatible with a granitic floor deepening towards the North in both cross sections. Over the La Maladeta granite the values of the anomaly are lower along profile AG1 than along profile MAL (Fig. 6). In cross-section MAL the observed gravity data above the La Maladeta granite are compatible with a granitic body 1 km-thick at the southern border that reaches almost 3 km depth in the North. In cross section AG1, the granitic body is larger, 2 km to almost 8 km-thick from South to North. In the northern end of cross section AG1 a minimum (−15 mGal) coincide with the outcrop of the Arties granite (Fig. 5b). The amplitude of this gravity low is compatible with a 6.5 km-deep elongated geometry for this granite as the calculated gravity data show after performing the 2.5D gravity modelling. In cross-section MAL, despite the absence of granite outcrops there, the observed gravity values are even lower (below −16 mGal) (Fig. 6). This gravity response is compatible with the presence of large volumes of granitic sills at depth as we have modelled to fit the calculated with the observed gravity data. This interpretation is based on the presence of the Bossòts Late Variscan granite (López-Sánchez et al. 2019) further north of the study area where big number of pegmatitic sills and veins crop out around of the main granitic body. In both cross sections, but more pronounced in cross section MAL, a relative small maximum, but still with negative values, appears at their northern sectors (Figs. 5,6). This feature has been interpreted to be associated with higher values of density of Devonian rocks affected by a higher metamorphic grade (staurolite-cordierite-andalucite zone) in that area (García-Sansegundo 1992;García-Sansegundo and Ramírez Merino 2013). The observed gravity data display an abrupt decrease coinciding with the beginning of the Nogueres Zone in both cross sections. This decrease reflects a wide negative anomaly of a total amplitude of ca. 15 mGal also oriented N120E (Fig. 3b) coinciding with the Nogueres Zone, the Ribagorza basin and the northern part of the Bóixols unit. It also correlates with the leading edge of the Rialp basement thrust sheet. In cross section AG1, the lowest gravity values (−18 mGal) coincide with the Aulet diapir (i.e. Triassic evaporites) (Figs. 2, 5b), whereas, in cross-section MAL, the lowest gravity values (also −18 mGal) do not correspond with Triassic evaporites outcrops at surface (Fig. 5a). The amplitude of this gravity low is compatible with a thick accumulation of Triasssic evaporites in the subsurface. In the southern end of the AG1 profile (Fig. 5b) there is a relative gravity maximum related to a duplication of Jurassic, Lower Cretaceous and Upper Cretaceous rocks due to the Boixols alpine thrust. At the southern part of both cross sections, the depth to the autochthonous deposits (aprox. 5 km depth) is based on previous geological cross sections by Teixell and Muñoz (2000) and Muñoz et al. (2018).

Discussion
Cross sections MAL and AG1 show good geological and geophysical constraints at their southern part. Moreover, some points have been refined through the gravity modelling, i.e. the distribution and thickness of the Triassic evaporites (poorly defined in seismic sections). The southern half of the Axial Zone, without available geophysical constraints such as seismic data or wells, has been modelled to infer the volume and geometry at depth of silicic granitic bodies whose density is generally lower than the host rocks (e.g. Vigneresse 1990) and testing the solution of gravity in the tectonic context of the Pyrenean Axial Zone (i.e. system of basement thrust sheets with a complex structure). Even though the density contrast between the granites and the basement is low (2.65 vs 2.70 g/cm 3 ) and more constraints are needed to reduce the uncertainty in their geometry, in absence of more geophysical data, gravity data have allowed validating our model of the geometry at depth of granitic bodies and to investigate along-strike structural variations inside the same basement thrust sheets.

Southern axial zone and the inferred geometry at depth of granites
The residual gravity anomaly along both cross sections highlights a different gravimetric response of the Orri basement thrust sheet with respect to both the Nogueres Zone and the thrust sheet outcropping to the North of the Gavarnie thrust. The outcrop of the Orri basement thrust sheet coincides with a gravity high that can be explained by the absence of the less dense Mesozoic rocks that are covering the Paleozoic basement in the Rialp and Ribagorza basement thrust sheets to the south (Figs. 5, 7). This assumption is compatible with the geometry of the Axial Zone proposed by the ECORS-Pyrenees profile dominated by an antiformal stack (Muñoz 1992).
Gravity data also reflect the presence of the La Maladeta and Arties granites and also suggest the presence of granitic rocks at depth (under Devonian and Silurian layers in the Aran Valley Syncline) related to the Bossòst granite located to the North (Fig. 7). With respect to the La Maladeta granite, both cross sections show a quite different geometry, the volume of the Maladeta granite in cross section MAL being considerably smaller than in cross-section AG1. The gravity values are compatible with a roughly laccolith geometry with a flat floor tilted towards the North in the cross section MAL, whereas in the cross section AG1 this granite is thicker with a steep floor tilted towards the North (i.e. towards the Gavarnie thrust). These differences could be related to the polydiapirism nature of the Maladeta complex (Leblanc et al. 1994) or simply to a possible originally irregular laccolith geometry. In addition, the cross section MAL crosses the Aneto unit, a well differentiated and concentric normally zoned unit, of the La Maladeta granite that could correspond to the external zones of the possible original laccolith. While the cross section AG1 crosses the larger and more complex zone of the granite, with no regular zoning and constituted by several coalescencing units (Leblanc et al. 1994). This zone could be interpreted as the central part of the laccolith. The deeper root of this granite, close to the Gavarnie thrust, could be associated with the interpretation of this thrust as an old major structure active during late Variscan extensional phases (Soler et al. 1998;Clariana et al. 2009;Teixell et al. 2018). During this period, large structures oriented in a roughly E-W direction have been interpreted as structural anisotropies that could facilitate the ascent of magma (Denèle et al. 2008;Clariana 2017). The laccolith geometry interpreted for the La Maladeta granite in cross section MAL could also be favored by the presence of Silurian shales, an effective décollement level introducing a rheological contrast that could help the lateral spread of magma. This behavior has been recognized in the neighboring Marimanha granite (Antolín-Tomás et al. 2009).

Evaporitic accumulation and underthrusting of basement thrust sheets
The residual anomaly map and the 2.5D gravity modelling of cross sections MAL and AG1 show a wide negative low (−17 and −18 mGal in cross sections AG1 and MAL, respectively) (Fig. 6). This gravity minimum shows an elongated geometry oriented N120E in the residual anomaly map 1 3 and parallel to the main structural trend (Fig. 7). It coincides at surface with large outcrops of Triassic evaporites belonging to the Ribagorza basin and at depth, with the hangingwall cutoff and leading edge of the Rialp basement thrust sheet, also oriented in an approximately N120E direction (Fig. 7). The amplitude of this gravity low is compatible with an important accumulation of Triassic evaporites at depth below the Jurassic and Cretaceous rocks in the Ribagorza basin and also in the Bóixols thrust sheet (Fig. 5).
Previous interpretations based on seismic data and wells have already considered an important accumulation of Triassic evaporites at depth in that sector (Muñoz 1992;Teixell and Muñoz 2000;García-Senz 2002;Mencos et al. 2015;Muñoz et al. 2018). Triassic evaporitic accumulations at surface have been interpreted as related with a centrifugal migration of the décollement level from the center of the South Pyrenean Central Unit towards its margins (squeezing) caused by differential loading and shortening during the Mesozoic-Cenozoic syntectonic filling (Soto et al. 2002;Santolaria et al. 2016). Cross sections MAL and AG1 show that the thickest evaporitic accumulation at depth is located above and in front of the leading edge of the Rialp thrust sheet. This feature suggests that the emplacement of the Rialp thrust sheet that occurred during the Oligocene (Muñoz 2002), could favor the accumulation of the décollement level at its southern leading edge. On the other hand, this area has been recently interpreted as part of a large diapiric province located along the Iberian Margin of the Pyrenean rift and defined by principal E-W and secondary N-S trending diapirs (Saura et al. 2016). These preferential orientations are also reflected by the pattern of the negative residual anomaly distribution over the northernmost part of the South-Pyrenean Zone studied in this work (Fig. 7).

Conclusions
This work shows that the previous proposed overall structural architecture of the Central Pyrenees is consistent with the observed gravity data in the study area, although some modifications (e. g. the shape and thickness of the evaporitic material) had to be introduced in order to better fit the gravimetric anomalies. The two 2.5D modelled cross sections show similar pattern in the observed gravity data that match well with similar interpreted structures in both cross sections. The observed gravity curves show a gravity high oriented N120E coinciding with the Orri basement thrust sheet and an important gravity depression, also oriented N120E, coinciding with the leading edge at depth of the Rialp basement thrust, which is interpreted in terms of accumulation of salt at base of the Mesozoic sequence. Another result of this work is the geometry at depth of the La Maladeta and Arties granites that have been constrained through gravity modelling. More geophysical data are needed to better constrain that geometry; nonetheless it is a first step towards characterizing the granites of the Central Pyrenees. Moreover, the feedback between geological and gravity data has been useful for constraining the thickness of the different Paleozoic units (Silurian, Devonian and Carboniferous) that form the basement thrust sheets in this part of the Pyrenees.
Acknowledgements This work is part of the project CGL2017-84901-C2-2-P funded by MCIN/AEI/10.13039/501100011033 and "ERDF A way of making Europe" and project PID2020-114273GB-C22 funded by MCIN/AEI/10.13039/501100011033 from Spanish Ministry of Science, Innovation and Universities. Seequent has provided us the GM-SYS module of the Oasis Montaj ® . The authors acknowledge the contribution of José María Llorente and Agustin González for the acquisition of the gravity data. We thank to Aigüestortes National park and Alt Pirineu Natural park their logistic support. We thank anonymous reviewer for improving the content in the manuscript. This study represents a contribution to GeoAp Research Group (E01-20R) (Aragón Government).
Funding Funding for this work comes from project CGL2017-84901-C2-2-P funded by MCIN/AEI/10.13039/501100011033 and "ERDF A way of making Europe" and project PID2020-114273GB-C22 funded by MCIN/AEI/10.13039/501100011033. Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.