3D Constrained Gravity Inversion and TEM, Seismic Reflection and Drill-Hole Analysis for New Target Generation in the Neves–Corvo VMS Mine Region, Iberian Pyrite Belt

Located in the Iberian pyrite belt, the Neves–Corvo mine is a world-class massive sulfide deposit and the largest operating mine in Portugal with underground mining down to 1000 m depth focused on massive and stockwork Cu, Zn, Pb rich ores. Gravimetric data have had a leading role in the discovery of the seven known deposits, together with time-domain electromagnetic (TEM) ground data. In this work, we present the results of a 3D constrained gravity inversion carried out with legacy ground gravity data. The 3D gravity inversions were carried out using an updated density database containing approximately 142,000 measurements. A recently constructed 3D geological model based on reprocessed 2D seismic reflection, 3D seismic, TEM and updated geology from detailed surface mapping and drill-hole data, was used to constrain the inversions. The results show multiple high-density anomalies that may indicate the presence of mineralization at depth. These anomalies were therefore cross-checked with holes previously drilled. Approximately 97% of more than 1000 available surface drill-holes located on or at a distance of less than 200 m from the high-density anomalies intersected mineralization. However, gravity anomalies have been drilled in the past and particularly dense black shales or rhyolitic/gabbroic rocks have been intersected. To increase the success of future drilling, gravimetric anomalies have been correlated spatially with high-conductivity TEM zones and strong-amplitude seismic reflections, because igneous rocks usually present weak-to-moderate conductivity and a massive column of black shales presents a seismic signature quite different from that of mineralization. We concluded that some of these locations represent high-quality targets to consider following up with drilling and further exploration.

Located in the Iberian pyrite belt, the Neves-Corvo mine is a world-class massive sulfide deposit and the largest operating mine in Portugal with underground mining down to 1000 m depth focused on massive and stockwork Cu, Zn, Pb rich ores. Gravimetric data have had a leading role in the discovery of the seven known deposits, together with time-domain electromagnetic (TEM) ground data. In this work, we present the results of a 3D constrained gravity inversion carried out with legacy ground gravity data. The 3D gravity inversions were carried out using an updated density database containing approximately 142,000 measurements. A recently constructed 3D geological model based on reprocessed 2D seismic reflection, 3D seismic, TEM and updated geology from detailed surface mapping and drillhole data, was used to constrain the inversions. The results show multiple high-density anomalies that may indicate the presence of mineralization at depth. These anomalies were therefore cross-checked with holes previously drilled. Approximately 97% of more than 1000 available surface drill-holes located on or at a distance of less than 200 m from the highdensity anomalies intersected mineralization. However, gravity anomalies have been drilled in the past and particularly dense black shales or rhyolitic/gabbroic rocks have been intersected. To increase the success of future drilling, gravimetric anomalies have been correlated spatially with high-conductivity TEM zones and strong-amplitude seismic reflections, because igneous rocks usually present weak-to-moderate conductivity and a massive column of black shales presents a seismic signature quite different from that of mineralization. We concluded that some of these locations represent high-quality targets to consider following up with drilling and further exploration.

INTRODUCTION
Several buried massive sulfide deposits in the Iberian Pyrite Belt (IPB), associated with the Volcano-Sedimentary Complex (VSC) (Fig. 1a) have been discovered in SW Iberia using the integrated interpretation of geological and geophysical data . This is the case of Neves-Corvo (Albouy et al., 1981;Leca et al., 1983), located south of Castro-Verde, in southern Portugal. Neves-Corvo is a world-class deposit with estimated massive sulfide tonnage totaling over 300 Mt. The seven known volcanogenic massive sulfide (VMS) deposits of Neves-Corvo are distributed along a NW-SE trend (Fig. 1b), dispersed in a large complex antiform structure (e.g., Carvalho et al., 1996;Oliveira et al., 2013) and associated with a large thrust fault, the so-called Neves-Corvo main (NCM) thrust (Pereira et al., 2021, Fig. 1a).
Understanding the subsurface deposit geometry requires an accurate geological map/model based on surface surveys, drill-hole logging and geological section interpretation . The accuracy of the conceived geological and ore deposit model at Neves-Corvo is however limited by the sparse geologic outcrops of host rocks. The model is also biased based on the distribution of drill holes. This opens the door for geophysical data to improve the geological modeling. Geophysical data acquisition, processing and interpretation allow costs reduction and a more regular sampling than drillhole data.
Gravimetry has historically been the geophysical method of choice in the region, due to the large contrast between the high-density massive sulfide deposits and the host rock lithologies (d > 4 g/cm 3 in massive sulfides and 2.9 < d < 3.9 g/cm 3 in stockwork veins, see Table 1). Host rocks include felsic volcanic rocks, black shales of the VSC (Neves Fm., see Geological Setting section) and shales and quartzites of the PQ Formation (Fm.), commonly presenting densities below 2.9 g/cm 3 (without mineralization).
Constrained gravity inversion has been performed in the past at Neves-Corvo for mapping new exploration targets (e.g., Lombador and Semblana; Lundin Mining Corporation website, Feb. 2016, h ttp://www.lundinmining.com/), though the results have remained unpublished. The construction of an accurate 3D structural geological model is essential to provide reliable targets. Incorporating new geological information and continually revising the 3D geological model make it worthwhile to revisit the constrained gravity inversion.
In this study, we carried out constrained 3D gravity inversion using legacy land gravity data, which have been revised and re-leveled. The constraints were imposed by a new 3D geological model based on reprocessed 2D seismic reflection profiles, time-domain electromagnetic (TEM) data, a 3D seismic volume, updated geological outcrop and cross section information (Matos et al., 2017, and up-to-date drill-hole information (Carvalho et al., 2020a. A recent density database with approximately 142,000 measurements (IPB LNEG density and Somincor/Lundin Mining Corporation databases; Marques et al., 2019) was used as the source for petrophysical constraint in the inversion.
The gravimetric inversion results are here investigated with the aid of TEM, seismic reflection, aeromagnetic and drill-hole data. Several high-density bodies produced by the inversion match highconductivity layers in the TEM 1D inverted cross sections. Some of these overlap with seismic reflectors. Several of the high-density anomalies were also intersected by drilling (or quite close to the drillholes) that have intersected mineralization. We believe these cases represent excellent targets for further exploration and drilling.

GEOLOGICAL SETTING
The study area is part of the IPB (Fig. 1a, b), which is considered to have been formed in an intraplate, continental, submarine volcanic environment during Paleozoic times (Oliveira et al., 2013;Inverno et al., 2015a). The Neves-Corvo region stratigraphy can be found, for example, in Oliveira et al. (2004Oliveira et al. ( , 2016Oliveira et al. ( , 2019 or Pereira et al. (2021). It comprises, from the oldest to the most recent, the following main geological formations (see Fig. 1c): (1) the Phyllite-Quartzite Group (PQ) basement, of middle Givetian to late Famennian age, with an unknown base (Mendes et al., 2020); (2) the VSC of Famennian-Late Visean age; and (3) the middlelate Visean Mé rtola Fm., a flysch sequence composed mostly of intercalations of greywackes and dark gray shales (Oliveira et al., 2004(Oliveira et al., , 2013, which is found in multiple thrusted sheets lying beneath the upper VSC. In this situation, the Mé rtola Fm. is usually designated as Mé rtola 2 and Mé rtola 3 (Oliveira et al., 2004;Pereira et al., 2021).
The VSC can be divided into two major sequences (Oliveira et al., 2013(Oliveira et al., , 2016. The upper VSC is composed, from bottom to top, of the Graça Fm. (black graphitic shales and gray siliceous shales) with mafic intrusive rocks and interbedded felsic volcanic rocks (Oliveira et al., 2013), the Grandaços Fm. (black shales with carbonate lenses and nodules), the ''Borra de Vinho'' Fm. (purple and green shales), the Godinho Fm. (volcanogenic sediments and gray siliceous shales) and the Brancanes Fm. (black shales pyritic shales and minor greywackes). The lower VSC autochthonous sequence includes the Corvo Fm. (black shales with tuffs and breccias), which rests in continuity on top of the PQ Basement, and the Neves Fm. (black pyritic shales), which is topped by a jasper and carbonate unit that represents the massive sulfide horizons hanging wall.
Two major VSC structures, the Rosá rio-Neves Corvo and Sã o Pedro das Cabeças Antiforms are present in the study area (Oliveira et al., 2013, 2016. These structures are located in a vast area where flysch sediments of the Mé rtola Fm. outcrop. The PQ outcrops in the Rosá rio-Neves-Corvo Antiform but it has not been detected so far in Sã o Pedro das Cabeças. The typical IPB southwest verging thrusts faults and nappes are also observed in the study area (e.g., Inverno et al., 2015a). The most important one, from a mineral exploration perspective is the NCM Thrust that breaks the two above mentioned major VSCs sequences (Fig. 1a) and tops the mineralization. The NCM thrust, its splays and other thrusts faults often cause the partial repetition/duplexing of the geological sequences, as previously mentioned for the Mé rtola Fm., but also of upper VSC sequence sediments and more rarely of Lower VSC sequence. The NCM thrust and the study area are all affected by Late Variscan nearvertical, SW-NE and N-S trending strike-slip faults, sometimes having a normal component (see Fig. 1b).

LAND GRAVITY AND DENSITY DATABASES
In the IPB Portuguese sector, the former government agencies Serviço de Fomento Mineiro (SFM) and Instituto Geoló gico e Mineiro (IGM), presently LNEG, have carried out systematic geophysical surveys, particularly gravimetric, in the Neves-Corvo region. These surveys were acquired during the second half of the twentieth century with the goal of locating exploration targets and attracting international investors. These surveys were complemented by others performed by Somincor, resulting in a reasonable gravity coverage across the IPB Portuguese sector (Represas et al., 2016;Matos et al., 2020).
Ground surveys over N-S, E-W grids, with distances between survey stations of 200, 100 and 50 m grid size were acquired between the 1960s and the 1990s (Oliveira et al., 1998). More than 388,000 gravity measurements have been collected in the IPB Portuguese sector. The Portuguese Geological Survey, LNEG, is responsible for the storage and maintenance of this database. The data were compiled and leveled using common data points. They were then gridded using a kriging algorithm using 500 m cells to achieve a uniformly distributed mesh. Smaller grid spacingÕs resulted in the generation of high-frequency signals with no apparent geological meaning. A regional Bouguer anomaly map of the IPB Portuguese sector has recently been produced (Represas et al., 2016) and interpreted ). The regional Bouguer anomaly was determined using a crustal density of 2.6 g/cm 3 , similarly to the density used in the IPB Spanish sector gravity surveys (Torres and Carvalho, 1998).
To ensure the quality and accuracy of the data, the legacy gravity dataset was reprocessed in this work to achieve a reliable 3D gravity inversion. To achieve this, it was essential to check the leveling of the different surveys and verify the consistency across the existing gravity database. For this purpose, a grid corresponding to the study area was extracted from the regional gravity map of the IPB Portuguese sector (Represas et al., 2016). The extracted leveled grid was compared to the original overlapping datasets to assess the presence of mathematical artifacts resulting from the leveling/gridding process. Some surveys, which were only available as paper maps, were digitized and the entire gravity surveys were re-leveled.
The resulting Bouguer anomaly map is shown in Figure 2a, with the data points of the different surveys overlain. Beyond the main anomaly, directly related to the Neves-Corvo VMS deposit (massive and stockwork ores), other geological features can be correlated with elements of the gravity map including VSC/PQ NW-SE lineaments and tectonic structures. Figure 2b shows the residual gravity map after the extraction of the regional field represented by a second-degree polynomial, with the geological contours shown in Figure 1b superimposed. Other methods of regional separation were used, such as other polynomials of different order and upward continuation to several altitudes but the second degree polynomial resulted in a better separation of the gravimetric anomalies. The residual map ( Fig. 2b) produced an enhanced separation of the gravimetric anomalies as it results in a signature that better preserves the signal assumed to be produced by near-surface geological structures of interest to this work.
The main thrust zones of the SW Neves-Corvo sector are represented by important gradients in the gravimetric maps. Subvertical, late, Variscan strikeslip faults of E-W, N-S and NE-SW orientations are also apparent from the residual anomaly map (e.g., N-S Lombador strike-slip fault). The relation of the gravity field with geological structures and massive ores is evident based on seismic reflection, drill-hole and geological data Matos et al., 2020;Carvalho et al., 2020Carvalho et al., , 2021.
To constrain our 3D inversion, the existing density database (LNEG and Somincor/Lundin Mining Corporation data) was updated to more than 142,000 density measurements. The Lower VSC units, poorly sampled in the past, have now been the focus of more than two hundred measurements. Table 1 presents rock densities measured from drill cores in the Neves-Corvo region. These values were used for the gravity inversions carried out in this study, averaged to the main geological formations used in the 3D geological model, i.e., Mé rtola Fm., VSC and PQ basement, plus the 7 known deposits of massive sulfides (see location in Fig. 1b).
Rock density is effected by several factors: porosity, fracturing, mineralogy, weathering and hydrothermal alteration. Considering the same geological formation, the presence of metallic minerals (e.g., sulfides or magnetite veins/disseminations) will increase its average density, while lowdensity minerals (e.g., clays, silica, sericite and carbonates) will have the opposite effect. The highest densities in samples with no apparent mineralization (3.22 g/cm 3 ) were found in the PQ basement, which can explain some positive gravity anomalies where the basement is shallower. These higher densities in the PQ basement may be due to the presence of hidden mineralization inside the core samples, as the PQ basement is sometime the source of stockwork mineralization and the vast majority of PQ samples present lower densities.

3D GEOLOGICAL INPUT STARTING MODEL AND CONSTRAINED GRAVIMETRIC INVERSION
The 3D geological model, used here as a parameter reference model in the gravimetric inversions, includes the mineralized deposits after drill-hole integration and the stratigraphic surfaces defined by the Mé rtola Fm.-VSC and the top of the PQ basement contacts. Both surfaces cover the entire study area, while the Upper and Lower VSC contact was also mapped, but only in the area of the known deposits, where a larger concentration of drill-holes with detailed geological logs exists  and, therefore, was not used in the present study.
After examination of the gravity database, average densities of 2.72, 2.77 and 2.79 g/cm 3 were established for the Mé rtola Fm., VSC, and PQ basement, respectively. Table 1 also shows us that Upper and Lower VSC present very similar average densities (2.76 and 2.77 g/cm 3 ), suggesting that the two units would not be distinguishable in the inversion process.
For the seven massive mineralization deposits, an average density of 4.53 g/cm 3 was used, while stockwork mineralization bodies were not considered in the parameter reference model. A look at the average densities shown in Table 1 shows that VSC host rocks and the PQ Fm. present similar densities and the greatest contrasts contributing to the gravity response are associated with the massive sulfides. The 3D gravity inversion was produced with a commercial software that uses a voxel-based approach (Seequent, 2020). The inversion code applies the Tikhonov minimum gradient regularization (Zhdanov, 2002) to produce a model that conforms to a given reference Earth model, and whose predicted response closely matches the input field data. The resultant model consists of a three-dimensional volume, discretized into a series of cells, each containing a predicted physical property.
The generalized inversion process can be represented by d = G m, where d is the data, represented in this inversion as the N residual Bouguer anomaly values, m is the model, corresponding to M density values and geological model interfaces, while G is the N 9 M inversion matrix.
The inversion process is carried out by iteratively minimizing the objective function uT = ud + k um, where ud is the data norm (difference between calculated and observed gravity field, designated as misfit), um is the model norm and k is the Tikhonov regularization parameter. k is set by the algorithm, and the regularization scale is made to vary in each iteration, so that the objective function is minimized and ud approaches 1.
The model norm can be expressed as um = W m (m 2 m0) 2 , where m0 is a reference model, and W m is a matrix of constraint weights, while the data misfit is given by ud = W d (d ÀGm 2 ), where W d are the a priori error levels defined by an uncertainty analysis and whose maximum values are set by the user.
The inversion code also makes use of the Cartesian cut cell method (Ellis and MacLeod, 2013) that allows a better representation of model interfaces, and of the iterative reweighting inversion focusing (Ellis, 2011), which is an iterative sharpening function that results in more confined anomalies that are geologically more plausible than those created by the smoothing constraints, necessary due to the under-determined inverse problem.
The voxel used as reference model, built from the 3D geological model ) have a horizontal cell size of 50 9 50 m and 25 m for the vertical, which is the geometry of the inversion voxel. A padding of 5 cells with a cell expansion ratio of 1.5 was added in all spatial directions. ure 3 displays the voxels used in the several inversion runs. The 3D geological reference model contains a high degree of realism, with the presence of several thrust faults and large folds (  Figure 1b and the location of the study area (black rectangle). A reference density of 2.6 g/cm 3 was used for the Bouguer anomaly calculation. Coordinates in Hayford-Gauss system, datum Lisboa, meters.
presence of the NCM thrust topping the Lower VSC and the mineralization) extends to the SE to the limits of the study area. To the NW, the Lower VSC and PQ Gr. are brought to the surface by the NCM thrust but the presence of stockwork mineralization in the PQ basement cannot be excluded. Table 2 summarizes the features of the three major runs carried out. In all runs, a weight constraint was set for each voxel of the reference model, representing the degree of confidence that the inversion should give to the respective value. This constrained some parts of the model, while allowing other parts to be calculated with a higher degree of freedom. Each of the three series of runs included two runs: (i) fixed limits and average densities for each layer and (ii) allowing the average density and layer limits to vary 20% (Table 1). The 100% constraint weight set for the first subrun is actually slightly changed by the inversion code, in order to accomplish the objective function criteria. In the second subrun, the constraint weight of 80% was chosen to allow density variations within the limits of the density database shown in Table 1. Several weight and smoothness parameters were tested for all the runs, but only minor variation of the density distribution were observed, probably because of the high geological constraints that were considered.
Despite the high level of constraints that were applied to the inversion, the absolute misfit (absolute difference between the calculated response of the model and the observed gravity field) was less than 0.04 mGal for all 6 runs, which means a relative error of less than 13%. A minimal difference in the inversion results and a misfit difference of 0.01 ± 0.005 mGal was observed between each of the subruns, that is, when a constraint weight of 100% or 80% were used. As expected, the lower misfits were achieved in the less constrained runs, with a minimum of 0.024 mGal (relative error of 8%) and these are the results presented here. Given the size and geological complexity of the inversion area, and the level of constraints imposed to the inversions, we found that these values represent a very good fit.
The result of the first series of runs ( Fig. 4a), which placed into evidence the importance of the massive sulfide lenses in the gravity field, considers d > 2.9 g/cm 3 . The inverted density model resembled the geometry of the known lenses, which was expected as these were fixed by 100% and 80%, but it should be noted that the inversion did not produce other high-density (above 3.1 g/cm 3 ) bodies. Densities between 2.7 and 3.1 g/cm 3 were placed surrounding the known deposits in both vertical and horizontal directions with thickness of 50-80 m, all below the strong Bouguer anomaly region. Also, in the area of high Bouguer values, the density anomalies located west of the known deposits were deeper than the known deposits. The northernmost one extended from 200 to 300 m to about 1400 m, while the southernmost one was centered at about 1.2 km, a bit deeper than the present-day exploitation depth of the Lombador deposit. The latter anomaly had a thickness around 300 m. Without any reference other than the deposits, the inversion process simply distributed densities < 2.7 g/cm 3 in a smooth way, filling in most of the inversion volume with densities around the reference density (2.6- Layer limits and average densities fixed 0.12 g/cm 3 initial contrast for Mé rtola with 100% constraint weight; 0.16 g/cm 3 initial contrast for VSC with 0.0001% constraint; 0.19 g/cm 3 initial contrast for the PQ with 100% constraint; 1.8 g/cm 3 initial contrast for mineral deposits and 100% constraint weight (i.e., a much larger degree of freedom is given to the VSC layer, where mineralized masses are known to reside, but also considering a high constraint weight for the already known deposits) Average densities and layer limits varying 20% 0.12 g/cm 3 initial contrast for Mé rtola with 80% constraint weight; 0.16 g/cm 3 initial contrast for VSC with 0.0001% constraint; 0.19 g/cm 3 initial contrast for the PQ with 80% constraint; 1.8 g/cm 3 initial contrast for mineral deposits and 80% constraint weight Density contrasts consider a reference density of 2.6 g/cm 3 2.7 g/cm 3 ), which was a relatively low density given the average density (not considering mineralization) in the region, according to Table 1 (2.76 g/cm 3 ). This run suggests that more high-density anomalies than the known deposits are needed to explain the observed gravity field.
In the second series of runs, the stratigraphic/ structural surfaces were spatially fixed and the VSC density was left unconstrained in order to see if the inversion converged to some of the known deposits. The results are shown in Figure 4b for d > 3.0 g/ cm 3 , and the anomalous volumes above the referred density have been marked with numbers and letters, for discussion purposes in the text. The inversion process has filled in the voxel's upper part (corresponding roughly to the Mé rtola Fm. and shallower part of the area of VSC outcrop) with densities between 2.55 and 2.65 g/cm 3 . The PQ basement was filled in with densities up to 2.80 g/cm 3 , as well as part of the VSC. Most of the voxel space was filled in with densities of 2.80 g/cm 3 . The VSC was also partially filled with densities that went up to the average density of massive sulfides (4.55 g/cm 3 ), although most of the VSC densities stayed below approximately 3.0 g/cm 3 , in a reasonable agreement with our database information for this unit (see Table 1).
In this inversion model, the densities ‡ 3.1 g/ cm 3 were spread across the study area but the largest mass concentrates were in the Neves, Corvo/Graça and Zambujal deposits area (marked with 1, 2 and 3, respectively, in Fig. 4b). However, these were pushed about 700 m south-westward of the true deposits position although still below the high Bouguer anomaly area (compare Figs. 2b and 4b). While most of the anomalies had thicknesses up to 400 m and depths above 1 km, this NW-SE oriented anomaly extended from 0.25 km down to 1.4 km, clearly below the deepest known deposit (Lombador deposit exploitation depth, which is opened at depth, is currently 1.1 km). This anomaly was placed over the NC thrust and is clearly associated with the difficulty of accurately representing the thrust in the reference model. Other reasonably sized high-density anomalies partially overlapped the Lombador and Semblana deposits (5 and 6, respectively, Fig. 4b), which is an encouraging result. One of the largest high-density anomalies was sited west of the Brancanes Anticline and seemed to have a structural imprint (marked 11, Fig. 4b), as well as another anomaly in the western limit of the study area (a, Fig. 4b), while other smaller anomalies were scattered across the study area. However, most of these smaller anomalies had densities between 3.0 and 3.2 g/cm 3 . Felsic rocks in the region have densities around these values. Anomalies with d > 3.2 g/cm 3 were limited to anomalies Testa (a), Gatoa (11), Neves, Corvo-Graça and Zambujal (1,2,3) and Caiada anomaly, located at the NE sector of the study area (d, Fig. 4b). Gatoa, Testa and Neves, Corvo/Graça and Zambujal anomalies had densities up to 3.7 g/cm 3 , which are well above that of local rocks without mineralization (see Table 1). In particular, Gatoa anomaly presented densities up to 4.3 g/cm 3 .
These results point to the importance of a good reference model to constrain the inversion, suggesting the importance of considering the known deposits. Furthermore, deviations of the modeled density anomalies from the true structures are expected as a result from the approximation of using a constant density for the Bouguer correction and the removal of a not truly real regional trend. The mineral deposits however, due to their high-density contrast with adjacent rocks, have an important impact in the gravity field and will contribute to deviations in the density solution, and thus should be included in the reference model.
The third and last series of inversions using the full 3D geological model, that is, the 3D structural model with the seven known deposits, (see Table 2) should therefore produce the most reliable results. The misfit of this third run was slightly higher than that of the previous runs, resulting from the higher level of constraint, with an average absolute misfit of 0.039 mGal. In this run, the density of VSC had a constraint weight of 0.0001% to allow the possibility of existence of other deposits. The inversion results are shown in Figure 4c for d > 3.5 g/cm 3 . It can be observed that the geometries of the known deposits were very well reproduced (they had close to 100% constraint) while several high-density anomalies that did not correspond to previously known deposits are also visible. These anomalies are marked with letters, which are related to their location. The anomalies detected in the second series of inversions (parameter reference model with only the 3D geological model) are attributed with the same letters here (a, 11, 1, 2, 3, d). Again, no important differences in the results were found when average densities and layer limits for Mé rtola Fm. and PQ Gr. were rigidly constrained or were allowed to vary by 20%.
It is important to note that, although the VSC layer, where the deposits are hosted, was left free, the lower densities in this layer output by the inversion result conformed to the lower limit of the density database ($ 2.5 g/cm 3 ). The densities below 2.6 g/cm 3 occupied a very reduced volume (approximately 5%) and were placed at the surface part of the VSC outcrop at the Rosá rio-Neves-Corvo and Brancanes Antiforms and some isolated places. These low density regions were spatially uncorrelated with the high-density anomalies (> 3.5 g/cm 3 ) discussed in the following paragraphs. The densities for the VSC layer varied between this value and that of massive sulfides (4.5 g/cm 3 ) but for the vast majority of the layer, the densities were between 2.74 and 2.8 g/cm 3 with minor regions and the upper part of the VSC outcrop presenting values between 2.6 and 2.74 g/cm 3 . These densities and the volumes they occupied were in agreement with the density database and the stratigraphic column. The Upper VSC had a large quantity of tuffs. This is the case of Godinho Fm., whose density measurements ranged 2.39-2.83 g/cm 3 in 477 samples in the database. Grandaços Fm. also bears volcanogenic sediments, with density ranging 2.58-3.01 g/cm 3 in 799 samples. In both formations, the upper limits were due to shales, while the lower bounds were due to tuffs. The nearly 2000 measurements in Mé rtola Fm. show densities varying 2.02-2.86 g/cm 3 . These values have good agreement with the densities output by the third run of inversions.
The vertical extent of the anomalies was generally 200-400 m for the isolated anomalies (thinner than the VSC thickness), while the central anomaly zone constituted by the Chanoca and Mestres-Algaré s anomalies had lower thickness, varying be-tween 50 and 200 m, in reasonable agreement with the thicknesses of the known deposits (i.e., 20-140 m). The lengths of the anomalies varied approximately between a few hundred meters (small Monte Branco deposits) and the 4-5 km of Chanoca and Mestres-Algaré s anomalies, that encompass a few of the known deposits. The depth of the anomalies varied between that of the known deposits (200 m) to approximately 1.5 km. The maximum depth of the known deposit is unknown, as the Lombador deposit dips to the NE and is open at depth, with current exploitation depth approaching 1.1 km.
The densities of the anomalies shown in Figure 4c ranged 3.5-4.5 g/cm 3 ; thus, in agreement with the densities of massive mineralization shown in Table 1. They were all placed inside the VSC layer, as this layerÕs density was set free and the other layers were not allowed to vary more than 20%; thus, could not host d > 3.5 g/cm 3 . In addition, we observed that the high-density anomalies, in particular those with d > 3.73 g/cm 3 , were placed close to the top of the lower VSC. This agrees with the known geological scenario for massive mineralization. Figure 5 shows the third series of inversion results for d > 3.7 g/cm 3 with the top of lower VSC surface overlaid. This surface was built using only drill-hole information because it was not possible to interpret this interface in the geophysical data. Therefore, the top of the lower VSC was only built in the central part of the study area where the deposits are sited and a significant number of drill holes is available. Nevertheless, this surface shows clearly that the high-density anomalies that may correspond to the presence of mineralization conform to the known geological scenario (massive + stockwork mineralization at the top of the lower VSC and near/below the main Neves-Corvo thrust).
We should also point out that most of the d > 3.5 g/cm 3 anomalies were placed below the region with high-Bouguer intensity (> 0.2 mGal) and analytical signal of the latter (see Fig. 10 in the Appendix). Though anomalies a, b, c, d, e, 8, 9 and 11 were slightly outside this region, they were still placed in an area with relatively high Bouguer anomaly (> À 0.1 mGal). These anomalies were also sited quite close to areas where the analytical signal is high (see Appendix), suggesting that they may be slightly misplaced due to errors in the reference model. This is particularly the case of anomaly b, which was sited over the NC thrust, which also caused the misplacement of the central anomaly in the second run.
To summarize, while the first run was capable of reproducing the observed gravity field with the presence of the known deposits and a few other high-density anomalies at the cost of a host volume with an average density value lower than the average of the measured density (see Table 1), the second run displaced the high density anomalies from the correspondent maximum Bouguer anomaly zone and their true position due to errors in the reference model and was incapable of replicating the geometries of the known deposits. Therefore, only a good reference model with the presence of the known deposits as in the third series of runs produced an inversion model with realistic densities for mineralization and local rocks, including that of host rocks (VSC), and compatible with the structure of the region, in spite of some errors in this reference model. In the following section, the high-density anomalies output by the third inversion run are investigated with the use of drill-hole, TEM and seismic reflection data, in particular the larger central area occupied by the Chanoca and Mestres-Algaré s anomalies. The first and second inversion runs were also checked against drill-hole and other geophysical data but the third run produced the best correlation with these data and for this reason it is the only one discussed below.

Investigation of Gravity Anomalies with Drill-Hole Data
Confirmation of whether or not high-density bodies (d > 3.5 g/cm3) may represent mineralization can be determined by checking the nearby drillholes. In areas where there is no such information, this could not be achieved. This was the case of anomalies Lombador North (9, Fig. 4c), Gatoa (11), and Testa (a) and their vicinities. The surroundings of Caiada anomaly (d) were investigated with four drill-holes, one of them very close to the anomaly, and no mineralization was intersected. The Cerro Grande anomaly (e, Fig. 4c) was also not investigated, although PCS-01 was drilled nearby.
Three high-density anomalies were intersected by drill-holes but no mineralization was detected. In particular, Castelo anomaly (b, Fig. 4c) was drilled by boreholes PSJ-43, PSG-33 and PSF-23-1 close to its southwestern limit, with no intersected mineralization. Anomaly Zambujal Sul (3b), which was crossed by drills SCA-14, SW-18, SR-20 and SU-20, all located close to the edge of the 3.5 g/cm 3 anomaly, also did not intersect mineralization. Anomaly Lombador NE (8 in Fig. 4c) was also crossed close to its edge by a barren drill-hole PNC-54. These three anomalies have densities compatible with the presence of stockwork mineralization but no mineralization was detected. It is possible that they were drilled too close to the edges of this body, within the error margins of the 3D input model. Less likely, but also possible, is that the non-intersection of mineralization might be due to drill-hole marker errors, as is the case for the drill-holes PSG-33 and PSF-23-1 that intersected anomaly Castelo (b), where the presence of stockwork and rubané (French term used to describe sheared stockwork subparallel to cleavage, thrusted over massive ore, Leca et al., 1983) was not considered in the lithologies database used in this work and was only detected after analysis of the original logs.
All the other high-density anomalies, including the ones not investigated by drill-holes mentioned above, remain potential targets. The top of the Castelo anomaly (b) was ''scratched'' by drill-hole PSG33 and intersected by drill-hole PSF33, where the presence of disseminated fine-grained pyrite in black shales and cherts and rubané had been detected, suggesting that stockwork mineralization may be the source of this anomaly. The anomaly disappeared for isoshell d > 4.3 g/cm 3 , supporting this interpretation. As admitted in the preceding section, this anomaly could also in part be due to the difficulty in representing accurately the NC main thrust in the reference model.
The d > 3.5 g/cm 3 anomaly of Chanoca (7, in Fig. 4c), that also occupied a large area between the Lombador and Semblana deposits and the area south and southeast of the latter deposit, was crossed by multiple drill-holes. West and south of Semblana deposit, several drills intersected massive mineralization, supporting the extension of this deposit in these directions (Pereira et al., 2021). Southeast of the Lombador deposit (5) and close to the Graça-Corvo deposit (2), several drill-holes crossed massive and stockwork mineralization. This suggests that at least a good portion of the Chanoca anomaly was caused by the presence of mineralization.
The Mestres-Algaré South anomalies were crossed by several drill-holes, of which seven intersected mineralization above the d > 3.5 g/cm3 anomaly at distances that varied from a few tens of meters to around 150 m. The Mestres-Algaré North anomaly was sampled above it and drill-hole PMM-02 detected mineralization 200 m above the anomaly, while drill PCA09-001, located 250 m north of the anomaly also crossed several levels of massive and stockwork mineralization. Therefore, we concluded that the large, central density anomaly d > 3.5 g/cm 3 , constituted by the Chanoca and Mestres-Algaré anomalies and that was placed below the region with high Bouguer anomaly and analytical signal (see Fig. 10 in the Appendix), generally represents the presence of mineralization and that further investigations are required to estimate the amount.
Finally, the Graça-Corvo South anomaly (2b) was traversed by drill-holes S0-03, SM-04 close to its edge and drill-hole SO-04 in the central part of the d > 3.5 g/cm 3 anomaly, although not reaching its d > 4.3 g/cm 3 nucleus. A few tens of meters above the d > 3.5 g/cm 3 anomaly and before its end, drillhole SO-04 crossed several meters of massive mineralization; thus, suggesting the presence of mineralization beyond the depth of the drill-hole and explaining the gravimetric anomaly. The western side of the Neves deposit was also investigated by multiple drill-holes and several of them had detected mineralization a few tens of meters above the d > 3.5 g/cm 3 anomaly.
To summarize, only 19 of the 661 surface drillholes that intersected mineralization were at a distance greater than 200 m of the d > 4.3 g/cm 3 gravimetric anomalies. Most of the available surface drill-holes were located over the known deposit and had obviously intersected mineralization but it should be noted that the known deposit geometry had a 20% degree of freedom in the inversion process and still they had a near 100% match with the geometry of the high-density anomalies output by the inversion. It should also be noted that mineralization was also detected in some drill-holes where no gravimetric anomalies exist. Further, five barren drill-holes crossed the high-density anomalies. Both situations are discussed in the section Integrated Analysis and Discussion below.
Of particular interest is the large central area occupied by the Chanoca and Mestres-Algaré d > 3.5 g/cm 3 density anomalies, which may be seen as an overestimation of the masses. However, mul-tiple drill-holes had crossed mineralization around the known deposits and disseminated pyrite was found everywhere in this region, suggesting that the results, even if somewhat overestimated, relate to the known geological setting. In the following two sections, the gravity anomalies are investigated with TEM and seismic reflection data in order to assess if they correspond to mineralization.

Comparison of High-density Anomalies with TEM Data
The high-density anomalies were next checked against TEM data, which were composed of 2D cross sections generated from the 1D inversion of ground-loop data with stations spaced every 100 m along each line and stitched together side by side. The data were acquired with a Crone PEM system between 2004 and 2012, with a transmitter loop size of 1 km 9 1 km and an average of approximately 13 Amps of current intensity, with 17 time windows of 0. 42, 0.56, 0.744, 0.988, 1.312, 1.744, 2.318, 3.078 4.088 5.428 7.206, 9.570, 12.66, 16.05, 22.7, 32.7 and 42.70 ms. The data were originally inverted using UBC EM1D code and several of the lines crossing the known deposits and high-density anomalies produced in the gravimetric inversion were reprocessed in the scope of this work. Beowulf software (Raiche et al., 2007) and commercial Maxwell software were used to invert the vertical component of the secondary field data.
In the Lombador area, the 1D inverted TEM data were constrained with about 300 drill-holes  but in the other areas where gravimetric anomalies were placed, drill-holes were sparse. Therefore, to perform the 1D inversion in this work, we used the following strategy. Two distinct and independent starting models were used to invert the data in the Lombador area, where multiple drill-holes are available; next, we checked which of these two initial models provided the best results when compared with the inversion constrained by drill-holes. The first starting model used was a 8layer model with 100 m thick layers over a halfspace, all with a constant resistivity of 100 X.m. The second starting model also comprised 8 layers over a half-space, each equally 100 m thick but with variable resistivities. The first layer and the half-space had initially the average resistivities found in the region for the Mé rtola Fm. and PQ basement, respectively, while the other 7 layers were attributed initial resistivities varying between 100 and 200 X.m. The two starting models were also used as reference models.
We verified that the second initial model produced the best results when compared with the constrained inversion and, therefore, this was the reference model used in the areas with little drillhole control where some of the high-densities anomalies were sited. The inversions were more dependent on the low resistivity layers present in the data (shallow aquifers, deposits or graphitic and pyritic shales), rather than on the model parametrizations or starting models. After the inversion, the 1D data points were interpolated/ smoothed between stations using EMITÕs Maxwell code to produce the 2D TEM sections. The interpolation process was carefully checked to avoid introducing false anomalies.
As stated above, several TEM lines that crossed the high-density anomalies produced by the 3D inversion were selected to be reprocessed using the process defined above: lines Semblana 2000S and 2600S, Guedelhinhas 2400 N, that intersected respectively, anomalies Zambujal Sul (3b in Fig. 4), Gatoa (11) and Mestres-Algaré North (10), and TEM lines Lombador 1100 N, 1700 N that totally crossed anomalies Lombador North (9), and lines Lombador 100 N and Lombador South 200, which intersected the Lombador deposit (5) and anomaly Lombador Northeast (8). The reprocessed TEM lines are shown in Figure 6, together with the density model. Their location can be seen in Figure 4c. The known mineralized deposits produced clear strongly conductive anomalies (white arrows, Fig. 6). This is clear in the original and reprocessed TEM lines although there was a better vertical coincidence between deposits/EM anomalies in the original 1D constrained drill-hole constrained inversion (not shown here). Although less intense as the mineral deposit, the TEM cross sections often presented several conductive layers.
An upper conductive layer (dashed arrows in Fig. 6) was often observed in the TEM sections, caused by water conducting thrust faults (fissural aquifers) or porous volcanogenic sediments in the upper VSC, or even the Brancanes Fm. carbonaceous, graphitic black shales with disseminated pyrite, mostly diagenetic Matos et al., 2020b). A regional conductive layer was also observed below the Neves-Corvo thrust fault (arrows in Fig. 6) and it is very possibly associated with the Neves Fm. black shales and local mineralization (e.g., the case of the Lombador deposit). This geological unit is one of the key ore horizons in the IPB Portuguese sector (Matos et al., 2011;Pereira et al., 2021). A near-match between the areas with highest conductivity in this regional conductive layer and the high-density anomalies (above 3.3-3.5 g/cm 3 ) can suggest the presence of mineralization. However, black shales in the IPB sometime present high- densities that may reach 2.87 g/cm 3 and, in some places, are the origin of gravimetric anomalies related with a significant volume of sediment (unpublished Empresa Mineira da Serra do Cercal report, 1996), or felsic volcanic rocks located in the lower VSC (West & Penney, 2017;Pereira et al., 2021). Igneous rocks generally present weak-moderate conductivity and can be distinguished from mineral deposit in TEM data. Black shales can, on the other hand, produce conductive zones as explained above, which at depths over 1.2 km where the EM signal is weak  can be difficult to separate from small-moderate sized mineral deposit signature. Figure 6 shows that the high value anomalies, with densities between 3.4 and 4.3 g/cm 3 , had good spatially coincidence with high conductivity regional EM anomalies (white arrows). The vertical position of these high-density/conductivity anomalies just above the PQ basement and within the top of Lower VSC represent a similar geological scenario to the Neves-Corvo mine known deposits, as stated in the previous section (see Fig. 5). The upper conductive zone generally did not present a high-density except in some areas where pyritic black shales were present; thus, the coincidence of a lower conductive zone with high-density is a good indicator of the possible presence of mineralization.

Gravity Anomalies and Seismic Reflection Data
The next step was to verify if the gravimetric high-density bodies were intersected by any of the available 2D seismic profiles or if they were sited inside the available 3D seismic volume (Yavuz et al., 2015). Massive mineralization and even stockwork mineralization have good contrast of acoustic impedance (AI) with adjacent rocks (Yavuz et al., 2015) and may provide strong reflectors if the mineralization thickness exceeds the vertical and lateral resolution of seismic data.
2D seismic surveys were carried out by Compagnie Gé neral de Gé ophysique (CGG) in 1991 and 1996, using explosives as the seismic source and 48 or 96 active channels. In total, 11 profiles were acquired (Carvalho et al., 1996). The 2D 1996 survey comprised six profiles that were recently reprocessed  with a clear improved imaging of the local structure. Profile 2 of the 1991 survey crossed one of the high-density anomalies, the Cerro Grande one (d > 3.4 g/cm 3 ), showing a set of linear seismic reflectors in the area occupied the gravimetric anomaly and suggesting that it might indicate the presence of the NCM thrust and associated mineralization. Profiles 4 and 6 of the 1996 survey intersected the Lombador and the Neves deposits, respectively (see Fig. 2, Donoso et al., 2020). The Neves deposit was crossed at the northeastern end of profile 6, where a weak to moderate reflector can be observed. The presence of the Lombador deposit is clearer in profile 4 , presenting a very strong reflector at its northeastern end, due to its large average thickness, which is around 100 m, attaining a maximum of 140 m in its central part. The central area of profile 6 traverses the Castelo anomaly in its northern part, coinciding with a strong reflector that presented a similar dip to that of the high-density anomaly. The central part of seismic profile 4 also intersected the Mestres-Algaré South anomaly, both the d > 3.4 g/cm 3 and the d > 4.3 g/cm 3 anomalies. A small but strong reflector coincided with the anomaly and multiple diffractions are visible, which may either come from faults or/and massive sulfides body.
One should bear in mind that depth conversion of the 2D seismic profiles was done with stacking velocities and velocity errors up to 20% may be present (e.g., Al-Chalabi 1979;Blackburn 1980), which result in depth positioning errors. Nevertheless, these encouraging results led us to inspect the available seismic volume (Ziramov et al., 2016), where depth conversion used not only stacking velocities but also drill-hole velocity log data from the Semblana region (Yavuz et al., 2015;West and Penney 2017). The recent reprocessing of the 1996 survey  and a recently acquired 3D seismic survey at Neves-Corvo over the Lombador deposit showed that mineralization could be mapped nicely, possibly some of the minor bodies that constitute that deposit Donoso et al., 2021). Figure 7 shows three depth slices from a migrated 3D seismic cube (see Fig. 4c for location). In 2011, 3D seismic data were acquired over an area of 25 km 2 covering known Neves-Corvo VMS deposits using mid-range vibrators (19 tons). Later on, in 2011 and 2012, three other 3D data sets covering areas of 42 km 2 and 10 km 2 were acquired and merged with the previous survey, creating a total of 87 km 2 within the study area. Details on the data acquisition and processing can be found in Yavuz et al. (2015), Ziramov et al. (2016) and West and Penney (2017). 3D Constrained Gravity Inversion and TEM, Seismic Reflection Figure 7. Correlation between 3D seismic cube and high-density anomalies obtained from 3D constrained inversion (runs# 3) for d > 3.5 g/cm 3 (orange with transparency) and d > 4.3 g/cm 3 (purple). Anomalies in green spatially match seismic reflector; other anomalies/deposits in black (see location in Fig. 1) The survey block that was made available to the authors and used in this work was the 2011 survey, which consisted of overlapping source and receiver grids that covered an area 6 km 9 4 km (NW-SE 9 NE-SW) with a bin width of 7.5 m 9 7.5 m. The processing did not preserve relative amplitudes (Yavuz et al., 2015;Ziramov et al., 2016), and dip moveout (DMO) and post-stack migration were applied to the data, contrarily to the full 3D survey, to which pre-stack time migration was applied (Ziramov et al., 2016). Using the frequency content in the data, stacking velocities obtained during processing and drill-hole velocity log data from the Semblana region (Yavuz et al. 2015;West and Penney 2017), we estimated a vertical resolution of 18 m according to the k/4 criterion for this survey, while horizontal resolution at depths of 0.5 km and 1 km were estimated as 135 m and 190 m, respectively. This 3D volume occupies all the known deposits central area, but unfortunately did not cover the area of the Testa, Gatoa, Monte Novo, Mestres-Algaré North, Castelo, and Cerro Grande anomalies (however, as explained above, the two last anomalies were covered by 2D seismic profiles). Depth slices at 750 m, 860 m and 1270 m below the surface are shown in Figure 7.
The upper depth slice (750 m) clearly shows the presence of the Lombador deposit and parts of the large d > 3.5 g/cm 3 Chanoca gravimetric anomaly. The Caiada anomaly was sited deeper but. interestingly, a strong reflection showed up very close to its location. The 860 m depth slice showed again the Lombador and now also the Semblana deposits. Parts of the Chanoca anomaly were also observed, while a negative amplitude reflection is now seen at the position of the Caiada anomaly, still located further deep. Finally, the depth slice at 1270 m did not show the presence of the Lombador and Semblana deposits anymore, which were located above, but a small part of the Chanoca anomaly can still be observed between the Lombador and Semblana deposits and northeast and southeast of Semblana. A strong positive and negative amplitude seismic anomaly is observed now at the d > 4.3 g/cm 3 Lombador NE high-density anomaly. The Caiada anomaly, present at this depth, did not produce any reflector. It is possible that seismic stacking velocities were inaccurate at this site and may indicate a wrong vertical positioning of the seismic reflector. For example, the seismic reflectors and the known deposits also had a vertical displacement. It also should be noted that the gravimetric inversion, even if constrained, does not provide an accurate depth positioning.
It should be noted that the shallower deposits of Neves, Corvo-Graça, Zambujal and Monte-Branco were not imaged in the 3D seismic volume. This can be due to the geometry used to acquire the seismic data that favored greater depths. However, the highdensity anomalies of Mestres-Algaré South, Zambujal South and Corvo-Graça South were placed at average depths between 700 and 900 m and this cannot be the cause for not producing clear seismic reflectors. Possibly their thickness was below seismic resolution, contrarily to what was suggested by the gravimetric inversion results or they did not correspond to mineralization.
The inline and crosslines of the 3D seismic cube were also inspected, as the depth slices may cause some difficulties to the correlation with the gravimetric and conductivity anomalies due to lack of an accurate depth conversion and gravimetric and EM inversion depth errors. The location of the inlines and crosslines discussed next can be seen in Figure 4c. It can be observed that below the Mestres-Algaré South anomaly, several strong reflectors are present, and the same applies to the Lombador North anomaly (pointed by arrows in Fig. 8a). The Lombador NE anomaly is also coincident with several reflectors at the depth level and below the anomaly (arrows in Fig. 8b).
Cross line 187 and neighboring crosslines correlated quite nicely with the Chanoca anomaly (Fig. 8c). Although the upper reflector may be caused by the NCM thrust (dashed small arrows), the reflector below suggests the presence of mineralization (arrows). However, this anomaly was drilled in multiple sites and only in some drill-holes mineralization was intersected, suggesting its origin may be in part due to particularly dense Neves Fm. black shales with high content of disseminated pyrite, as no felsic or mafic intrusions are inferred from the total magnetic intensity map (not shown here, see Marques et al., 2020b;Matos et al., 2020b). The small, d > 3.5 g/cm 3 Caiada anomaly corresponded in the crosslines with a moderate reflector.
The inlines imaged the Chanoca, Lombador NE and Lombador N anomalies quite nicely (Fig. 8c). The latter, similarly to the Caiada anomaly, seemed to be vertically displaced from the reflectors around inline 143, which might again be due to inaccurate depth-conversion seismic velocities. The multiple reflectors observed below the gravimetric anomalies might be due to several levels of mineralization or dense volcanic rocks (rhyolites). Felsic volcanic rocks have been identified as the source of some seismic reflectors (West and Penney 2017).

Integrated Analysis and Discussion
As stated at the end of Sect. 3D Geological Input Starting Model and Constrained Gravimetric Inversion above, in all of the three gravimetric inversion runs, the third run (reference model with deposits and structural model) results (> 3.5 g/cm 3 density anomalies) produced the best correlation with mineralization detected (or not) in drill-holes and spatial coincidence with high-conductivity anomalies in the TEM sections and seismic reflectors. It was also the run whose density distribution and range of values better matched the known stratigraphy and the density database with more than 142,000 measurements. The following discus- sion is therefore on the results of our preferred model, i.e., the third inversion run.
The third inversion run was also the one that produced the largest volume of high-density anomalies. In some cases, denser black shales may also cause the gravimetric anomalies that can show up for d > 3.5 g/cm 3 . Looking at the TEM lines, the black shales may be more conductive than volcanic rocks, particularly if they contain disseminated pyrite. If one observes the TEM cross sections at these locations ( Fig. 6 and other TEM lines not shown here), it can be concluded that some of the highdensity anomalies spatially coincide with a peak of low resistivity at or nearby the TEM cross sections.
The Gatoa anomaly produced a positive anomaly but without very low resistivity values in the original 1D inversion. In our reprocessing, it produced a very low resistivity anomaly spatially coincident with a gravimetric anomaly, but vertically displaced about 100 m (white arrow, Fig. 6a). This may be due to the great depth ($ 1350 m). According to the reprocessing done here and elsewhere , the EM signal below 1200 m was weak. The Mestres-Algaré North d > 4.3 g/cm 3 and d > 3.4 g/cm 3 anomalies, which had no seismic coverage, produced in the TEM section that traversed it a high conductivity anomaly with a similar area, although slightly displaced from the high-density anomaly. The d > 4.3 g/cm 3 anomaly of Mestres-Algaré South, which had a seismic signature in the 3D volume, had a coincident EM low resistivity peak.
The Chanoca anomaly, which in many parts was imaged by the 3D seismic volume, also corresponded in parts to areas of high conductivity, although a vertical mismatch of about 150 m was observed between the anomalies. This was the case of the area between the Lombador and Corvo-Graça deposits, where several drill-holes intersected mineralization at the levels of both anomalies. To the southeast, between the Corvo-Graça and Semblana deposits, several drill-holes crossed the Chanoca anomaly, but until this time no mineralization was crossed, except west and southeast of Semblana.
The regional conductive layer (Neves Fm.) accompanies the Chanoca d > 3.5 g/cm 3 high-density anomaly (with about 100-150 m vertical displacement) until Semblana deposit. In the original inversions, in terms of conductivity, there is no significant difference in the values observed in the area where drill-holes were productive or where drillholes were barren. However, in the reprocessing of TEM data carried out in the scope of this work, the Semblana deposit presented a clear very low resistivity signature, again with vertical mismatch of approximately 120 m (see Fig. 6b). Therefore, the reprocessing of all TEM cross sections in this area or its 3D modeling can provide important clues to discriminate between productive and less potential areas of high-densities anomalies.
Finally, we analyzed the Lombador N and Lombador NE gravimetric anomalies in terms of conductivity response in the TEM data. TEM lines Lombador north 1700 and 1100 crossed the Lombador N anomaly (Fig. 9, see location in Fig. 4c). The central part of the anomaly was sited at 1350 m depth and in the original inversion, a conductive region was observed about 250 m below the highdensity anomaly. The two lines reprocessed here (we recalled it is an unconstrained 1D inversion) show a clear high conductivity anomaly placed between 250 and 350 m above the high-density anomaly (black arrows, Fig. 9). Taking into consideration the previous vertical mismatch between gravimetric and EM anomalies, and also the positioning of the seismic reflectors observed at and below the gravimetric anomaly (arrows in Figs. 8a, c and 9), this location constituted a very attractive target for future drilling. The Lombador NE anomaly was sited slightly above the Lombador N anomaly, at approximately 1250 m. TEM lines Lombador North 100 and 1100 had intersected the anomaly, and both were reprocessed (see Fig. 6c). A very high-conductivity anomaly can be observed with the d > 3.5 g/cm 3 anomaly, and in one of the TEM lines a peak of the high-conductivity anomaly matched the d > 4.3 g/ cm 3 anomaly (black arrows in Fig. 9). The match with the seismic reflectors observable in the seismic cube is also noticeable (blue arrows in Fig. 9).
To conclude, the presence of several high-density anomalies deserve further investigation. Obviously, the anomalies that were crossed by nearby drill-holes and detected mineralization or the anomalies that coincided with high-conductivity anomalies and seismic reflectors have the highest interest.
However, this analysis will be affected by errors intrinsic to the processing of gravimetric, TEM and seismic methods. As stated above, neither the gravimetric nor the EM methods provided accurate depths, while the seismic method without detailed velocity information suffered from the same limitation. We must also note that we used 1D inverted TEM data to investigate the gravimetric anomalies in a very complex and 3D geological medium. In spite of the fact that recent 3D forward modeling carried out in the Lombador area by some of the authors  did not significantly change the 1D inversion results, we should expect that 1D inversion in such a complex medium and the limited lateral extent of the mineral deposits will result in position/conductivity errors.
In particular, the gravimetric method, which was the main focus of this work, was affected by the non-uniqueness problem. Although we partially overcame the problem with an updated 3D initial reference model and new LNEG and Somincor/ Lundin Mining Corporation rock density database, inaccuracies in the model undoubtedly remained that will affect our results. Furthermore, an underestimation of the PQ basement density might have occurred, in particular in deep unknown areas (so far not drilled). We used an average value of 2.79 g/ cm 3 derived from more than 400 samples collected from drill-holes (Table 1). As drill-holes usually did not go deep inside the PQ, the value used here might indeed be underestimated. Alternatively, a denser pre-PQ basement might be present in some regions in the lower part of the model and which was not accounted for. Some of the reflection seismic profiles supported this hypothesis (e.g., Carvalho et al., 2017, Neves-Corvo SE sector). In spite of the large number of drill-holes in the area of the 3D inversion, inaccuracies in the top of PQ model may also be present, which together with the density underesti-mation may produce a high-density VSC. These two problems may contribute to the large central, d > 3.5 g/cm 3 anomalous zone constituted by the Chanoca and Mestres-Algaré anomalies, which, despite the presence of mineralization and disseminated pyrite observed in drill-holes, may suffer from some overestimation.
Therefore, the exact geometry of these previously unknown central zone and localized/isolated high-density anomalies is not entirely trustworthy. Furthermore, because of the mesh vertical cell size used (25 m) and the typical thicknesses of known deposits (vary between 15 and 30 m) that requires a model with thinner mesh cell (voxel) size, the geometries of the density bodies should be treated with caution. Also, the spacing between gravity data points poses a limit to the resolution capability of the dataset. In spite of these limitations, some of the high-density anomalies remain viable targets for drilling and exploration, as they are placed at intermediate depths (700-1200 m) and the d > 4.3 g/cm 3 anomalies have dimensions between 0.3 and 0.6 km, similar to the smaller Monte Branco and Zambujal satellite deposits.
The success of the gravimetric inversion can be assessed preliminarily by comparing the spatial coincidence of the gravimetric anomalies with existent drill-holes that crossed mineralization. Using only surface drill-holes, as the underground drill-holes were only performed over the mineral deposits, we Figure 9. Correlation between 3D seismic cube, TEM reprocessed data (see location in Fig. 4c) and high-density anomalies of Lombador North obtained from the 3D constrained inversion (runs# 3) for d > 3.5 g/cm 3 (orange with transparency) and d > 4.3 g/cm 3 (purple). TEM lines Lombador North 1700 and 1100 are shown. Arrows point to seismic reflectors (blue) and high-conductivity anomalies (black) associated with Lombador N gravity anomaly.
verified that only 19 out of the 661 drill-holes that intersected mineralization were at a distance farther than 200 m of the d > 4.3 g/cm 3 gravimetric anomalies (hereafter always referring to the third run), meaning that 97% of drill-holes with mineralization were inside or very close to the gravimetric anomalies. This percentage included drill-holes that intersected the known mineral deposits and that were used to constrain the inversion, but, as shown in Table 2, the known mineral deposits geometry was allowed to vary by 20%. On the other hand, only 4 out of 7 drill-holes that crossed the d > 4.3 g/cm 3 gravimetric anomalies that did not correspond to known deposits were barren drill-holes. We further point out that some of the density anomalies not associated with known mineralization remained consistent between iterations utilizing reference models both with and without the massive sulfide bodies.
Most of the outliers (gravity anomalies with barren-drill-holes and drill-holes that intersected mineralization that did not cross gravity anomalies) can be explained by errors in the 3D geological model used to constrain the inversion and/or limitations of the inversion process (e.g., non-uniqueness of potentialfield data). The use of EM and seismic data helped us to understand these few observed outliers and better assess if gravity anomalies do correspond to mineralization. Drill-holes with mineralization that did not coincide with d > 4.3 g/cm 3 gravity anomalies, except south of Semblana, were associated in the 3D seismic cube with weak-moderate reflectors of limited areal extension, suggesting that we might be in the presence of small mineralized bodies, which cannot be detected by the gravimetric data used. These locations were not covered by TEM data but nearby lines suggest they spatially match the regional conductive Neves Fm., layer confirming the possibility of small mineralized bodies. On the other hand, the scarce (4) barren drillholes sited inside a d > 4.3 g/cm 3 gravity anomaly were also not covered by TEM data but, again, they seemed to be associated with the mineralized conductive Neves Fm. The 3D seismic volume did not present any coherent reflector at these locations but it should be emphasized that it also did not image the Corvo, Graça or other of the shallow deposits, possibly due to acquisition geometry/statics problems.
Mineralization of important dimensions always produces conductive anomalies and strong reflectors in seismic data if the thickness is greater than the vertical and lateral resolution and appropriately acquired/processed, as shown by recent seismic surveys acquired in Neves-Corvo Do-noso et al., 2021). In spite that other lithological/ structural features may also produce conductive zones and seismic reflectors, the simultaneous presence of electrical, seismic and gravity anomalies at coincident locations provides additional confidence to the gravimetric inversion results. Some interesting gravimetric anomalies drilled in the past have proven to be thick, particularly dense back shales columns or rhyolitic rocks. Some of these sites may also produce seismic reflectors or conductive areas but more rarely they produce both. Therefore, the fact that some of these gravity anomalies obtained in this work also corresponded to high-conductivity anomalies, seismic reflectors, and mineralized (or nearby mineralized) intersections, increases the confidence in the gravimetric inversion results in these cases, although in some cases the presence of graphite and disseminated pyrite in shales can also be the causative source of the anomalies. This combination of gravity, seismic and EM results has led in the past to the discovery of deeper-seated deposits, like Semblana (Pereira et al., 2021).

CONCLUSIONS
In this work, legacy land gravity data were used to carry out a constrained 3D gravimetric inversion. The gravity data were reprocessed and a new Bouguer anomaly map was produced. A new 3D geological model, in conjunction with an updated petrophysical reference database, and reprocessed seismic reflection (2D and 3D) and TEM data, were used to constrain inversion and guide our interpretation of the results.
Several high-density anomalies resulted from the 3D inversion procedures. Due to non-uniqueness of gravimetric data and because high-density anomalies can also be produced by large volumes of basic rocks or black shales, these anomalies were investigated with seismic reflection, TEM and drillhole data, once more. Rarely, those rocks produce simultaneously high-conductivity, gravity and seismic anomalies, and their combined interpretation can contribute to distinguishing between mineralization and other geological causes. It was verified that several of the high-density anomalies were spatially coincident, or are close to nearby seismic reflectors and conductivity anomalies. Ninety seven percent of surface drill-holes where mineralization was intersected were sited nearby (in or at a distance less than 200 m) the high-density anomalies, making these anomalies very interesting targets for exploration and for future investigation with drill-holes.
An iterative approach that incorporates new or reworked geologic and geophysical datasets, as constraints to the inversion process, is the best solution to the problem of non-uniqueness in gravity inversion. This is true not only for gravity, but for all other potential field geophysical datasets. Adopting this type of adaptive methodology will lead to better results and the potential to reveal small and/or deeper orebodies previously undetected. The use of a finer inversion mesh will allow a better definition of geometry, but this would also require a higher resolution gravity survey. Improvement of the knowledge of the 3D model used to constrain the inversion and, in particular, the deeper structure and densities, or a constrained inversion using drill-hole densities at the respective drill's locations can improve the inversion results in future.

AVAILABILITY OF DATA AND MATERIAL
Gravimetric data are commercially available at LNEG; density database is available upon request.

DECLARATIONS
Conflict of Interest None.