The implications of crustal architecture and transcrustal upflow zones on the metal endowment of a world-class mineral district

Earth’s mineral deposits show a non-uniform spatial distribution from the craton-scale, to the scale of individual mineral districts. Although this pattern of differential metal endowment is underpinned by lithospheric-scale processes the geological features that cause clustering of deposits remains enigmatic. The integration of geological and geophysical (seismic, gravity, and magnetotelluric) features has produced the first whole-of-crust image through an iconic Neoarchean volcanic complex and mineral district in the Abitibi Greenstone Belt, Superior Province, Canada. Observations indicate an asymmetry in surface geology, structure, and crustal architecture that defines deep transcrustal magmatic-hydrothermal upflow zones and the limits of the Noranda District ore system. Here, extreme volcanogenic massive sulfide (VMS) endowment is confined to a smaller area adjacent to an ancestral transcrustal structure interpreted to have localized and optimized magmatic and ore forming processes. Although lithospheric-scale evolutionary processes might act as the fundamental control on metal endowment, the new crustal reconstruction explains the clustering of deposits on both belt and district scales. The results highlight a strong magmatic control on metal and in particular Au endowment in VMS systems. Overprinting by clusters of ca. 30 Ma younger orogenic Au deposits suggest the ore systems accessed an upper lithospheric mantle enriched in Au and metals.


Supplementary Appendix DR2. GEOPHYSICAL METHODS
A regional, deep seismic reflection profile was designed to image the structural architecture of the crust and covered a ~50 km transect centered on the Noranda District (Fig. 2b). Cheraghi et al. (2018) and Naghizadeh et al. (2019) presents the specifications of the survey and the processing performed to produce the migrated profile in Figure 2. In summary, an array of four vibrator trucks producing a linear upsweep of 2-96 Hz was repeated four times at each nominal source location. Surveying used 50 m source (4 sweeps) and 25 m receiver intervals. Acquisition parameters were a 12 or 16 s recording length, 2 ms sample rate, 15 km-0-15 km spread size, 5 Hz geophones (single), and AHV-IV 364 Commander vibrator trucks . The seismic data processing stream included trace kills and reversals, minimum phase conversion, ensemble balance and amplitude recovery, surface consistent scaling, linear and erratic noise attenuation, surface-consistent deconvolution, anomalous frequency suppression, refraction statics, velocity analyses, surface consistent residual statics, and postand pre-stack time migrations . There is some overlap with seismic transects (line 21 and 21-1) conducted across the Superior craton ~30 years ago as part of the Lithoprobe program (Calvert & Ludden, 1999;Percival & West, 1994;White et al., 2003). Overlaps between the two surveys show better resolution of the upper crust along the Metal Earth seismic profiles .
Gravity data were acquired from transect stations with an average spacing of 300 m and side road stations up to 1 km from the transect (Maleki et al., 2020). The data were combined with existing Geological Survey of Canada regional compilation (http://gdrdap.agg.nrcan.gc.ca/) and are corrected for changes in latitude, elevation using software Oasis Montaj® (https://www.seequent.com/productssolutions/geosoft-oasis-montaj/). Terrain corrections were applied assuming background density of 2.67 g/cm 3 . The VPmg software Version 9.0 (Fullagar and Pears, 2007;Fullagar, 2013) linked to Mira Geoscience GOCAD mining suite was used to perform a smooth unconstrained 3D gravity inversion. The initial 3D density model was constructed within a volume with dimensions of ~(84 km x 104 km x 15 km) and a 450 m cell size. Topographic information (SRTM 90m) used for this inversion was downloaded from Geosoft public DAP server (http://dap.geosoft.com). The inversion started with initial density contrast value of 0.01 g/cm 3 for all the cells in the subsurface model and allowed VPmg to numerically change the property of each cell to fit the observed data. The lower and upper density contrast values are fixed during inversion i.e. 0 and 1. The final unconstrained inversion of gravity data, after 78 iterations, returned the RMS misfit of 1 mGal from an initial value of 10.68 mGal. The density cross-section ( Fig. 2c) was extracted along the transect from the final inverted 3D density model. A 3D isosurface model (Fig. S2) was generated using the density contrast value of 0.07 g/cm 3 . The density cross-section delineates several low-density anomalies (G1-G6) and corresponding 3-D isosurface model provides geometry of low-dense features. The smooth unconstrained potential-field inversions suffer with problem of non-uniqueness (Fullagar and Pears, 2007;Li and Oldenberg, 1998) and hence should be combined with other geophysical methods for a plausible geological interpretation.

Metal Earth contracted the collection of 750 new MT measurements through Complete
Magnetotelluric Solutions and Moombarriga Geoscience, comprised of a combination of broadband-and audio-MT soundings using the Phoenix Geophysics V5-2000 MT system. When the newly collected MT data set is combined with the existing MT data sets there are upwards of 1100 MT soundings in the area of interest. A total of 13 'regional' north-south MT transects ranging in length from 60-150 km were acquired (Hill et al. 2021;Roots et al., 2022), which when combined with existing data sets such as Lithoprobe (Jones et al., 2014) provide an areal coverage of > 2.5×105 km 2 . Nested within these 'regional' transects are more densely sampled high-resolution segments up to ~10 km in length to provide improved discrimination of lateral upper crustal structure in zones of localized strain that are commonly coincident with mineralised domains. Broadband-MT data spans the frequency range of ~320-0.001 Hz, depths of up to ~75 km given favourable survey conditions (Chave and Jones, 2012). While the audio-MT data spans the frequency range, as the name implies, from the audio band in which sound waves are detectable to the human ear (~10000-1 Hz), a depth range of ~0.5-3 km can be imaged given the resistive near surface (Chave and Jones, 2012). MT soundings were collected primarily along profiles; however, a number of soundings were collected offset from the main profile trace. These stations, combined with the pre-existing data, comprise a 'swath' or 'corridor' transect, giving improved 3D control. The regional transects have a station spacing of ~5 km and are comprised solely of broadband MT soundings with observation recordings of 24-48 hours. The high-resolution segments within the regional transects contain alternating broadband-and audio-MT soundings with a station spacing of ~330 m and recording times of ~24 hours for broadband MT soundings and ~4 hours for audio-MT soundings. All time series data collected within Metal Earth were processed to obtain estimates of the frequency domain transfer functions Z and K using robust remote reference processing (Jones et al., 1989)  information that allows identification of structure away from the main transect trace. This off-transect 3D structure is required to honour the 3D nature of the observed data identified by phase tensor analysis.
Longer wavelengths probe deeper and farther from the observation location; however, the resolution correspondingly lessens with both depth and distance (Chave and Jones, 2012). Inverse modelling of all of the 'swath' or 'corridor' transects in the Upper Abitibi region supplemented by extant Lithoprobe data (Jones et al., 2014) allows treatment of the meandering profiles including off-transect soundings. HexMT (Kordy et al., 2016a,b;Wannamaker et al., 2017) was used to compute the inverse models presented.
Inversion was conducted on 112 sites with 18 logarithmically spaced periods between 144-0.0011 Hz and error floors on the real and imaginary parts of the complex impedance elements Zij applied of max (assigned Zij error, 5%(|Zxy-Zyx|/2)) at each frequency. The starting and apriori model for inversions was a 600 Ωm half-space mesh consisting of 207x237x59 nodes and a nominal horizontal cell size of ~1.2 km with areas of local refinement. The preferred best fit model with static shift solution had an nRMS of 2.18. The displayed sections follow the meandering nature of the measured Rouyn-Noranda MT soundings. The 3D inversion closely reproduces all modelled components of the observed data.
The HexMT algorithm (Kordy et al., 2016a,b;Wannamaker et al., 2017), exploits deformable hexahedral finite elements (FE) to precisely incorporate topography. It uses direct solvers parallelized on symmetric multiprocessor (multi-core) workstations with large RAM for both the forward solution, parameter Jacobians, and the parameter model update. First-order edge elements are used to represent the secondary electric field (E), yielding accuracy O(h) for E and its curl (magnetic field). A divergence correction for the E-field is applied using the single-step Hodge decomposition. The system matrix factorization and source vector solutions are computed directly using the MKL PARDISO library. The factorized matrix is used to calculate the forward response as well as the Jacobians of electromagnetic field and MT responses using the reciprocity theorem. By exploiting the data-space approach, the computational cost of the model update is greatly reduced in both time and computer memory compared to that of the model-space forward simulation. In order to regularize the inversion using the L2 norm of the gradient, we factor the matrix related to the regularization term and apply its inverse to the Jacobian, which is also performed using the MKL PARDISO library. For dense matrix multiplication and factorization related to the Gauss-Newton model update, we use the PLASMA library. The libraries show good scalability across processor cores.     Caption: Figure S4. Geophysical data locations map. Yellow color represents MT sites, red and black colors represent gravity locations, and seismic profile. The Seequent software Geosoft Oasis Montaj (https://www.seequent.com/ products-solutions/geosoft-oasismontaj/) was used to generate the data location map.