Constraints on Growth and Stabilization of the Western Superior Craton From Inversion of Magnetotelluric Data

A data set consisting of 376 broadband and long‐period MT measurements was used to generate the first ever 3D resistivity model of the Archean western Superior Craton. The modeled resistivity structure is compared to coincident seismic reflection data. The observed geophysical signatures are interpreted within the context of the late stages of crustal growth and cratonization of the region via the progressive accretion of terranes against the initial cratonic core. The northern‐most terranes comprising the cratonic core exhibit a nearly homogenous highly resistive crust. The lower crust of the southern terranes contains largely continuous low resistivity bands which run subparallel to major terrane boundaries and corresponding fault systems. In some cases, low resistivity features are coincident with dense packages of sub‐horizontal to listric reflections within the mid‐ to lower crust. These resistivity structures are inferred to represent preserved geoelectric signatures of late to post‐orogenic magmatic pulses likely related to delamination of locally overthickened crust. Increased mantle heat flow resulted in partial melting of the lower crust and upper mantle and upward migration of CO2‐rich melts and fluids through crustal weak zones corresponding to shear and/or suture zones formed during terrane amalgamation. Thermal softening of the mid‐ to lower crust led to orogenic collapse and reactivation of the crustal shear zones, resulting in formation and interconnection of graphitic films which were preserved within the stable craton. These results have implications for the tectono‐magmatic history of the western Superior Craton, as well toward the understanding of the geodynamic regime of the Archean Earth.


Introduction
Archean cratons form the stable cores of the Earth's continents (Hoffman, 1988), and the processes which led to their formation and long-term stability remain enigmatic given the difficulty of obtaining evidence of these ancient processes in modern geophysical images and the lack of modern analogs.Shifting rheological conditions during the Archean likely led to a secular evolution in the geodynamic processes responsible for crustal formation and stabilization, including the intermittent predominance (e.g., Dhuime et al., 2012;Perchuk et al., 2020) and/or transition to modern-style subduction-accretion systems (e.g., Chowdhury et al., 2021;Mole et al., 2021;Van Kranendonk, 2010).Furthermore, forms of lithospheric recycling which are uncommon in modern orogens may have occurred more frequently during the hotter Archean geothermal regime (e.g., delamination, shallow slab breakoff; Chowdhury et al., 2021;Johnson et al., 2014).These geodynamic processes, and the associated magmatic-hydrothermal activity and lithospheric architecture in turn have implications for the volatile budget of the resulting crust (Bierlein et al., 2006;Groves et al., 2016).Thus, constraints on the processes active in the Archean are critical for our further understanding of the formation and stabilization of these terranes and the enrichment/depletion and localization of volatiles (including base and precious metals) within Archean lithosphere.
Geophysical methods are capable of imaging subsurface architecture across hundreds of kilometers, at depths ranging from the upper crust to the mantle.In the last decade, the magnetotelluric (MT) method has emerged as a particularly useful tool for investigating such large-scale lithospheric features.Continental-scale studies across Australia (AusLAMP, Robertson et al., 2016), China (Sinoprobe, Dong et al., 2016), the USA (USArray, Bedrosian & Feucht, 2014;Kelbert et al., 2019), and Canada (Lithoprobe, Jones et al., 2014;Metal Earth, Hill et al., 2021;Roots et al., 2022) have proven crucial in elucidating the lithospheric architecture of the respective continents and geodynamic processes by which they formed and evolved.MT data is sensitive to subsurface resistivity (or its inverse, conductivity), which depends on a variety of factors.Within the lithosphere, these include temperature (e.g., Heinson et al., 2018;Selway, 2014), increased iron or fluorine content (e.g., Evans et al., 2011;Robertson et al., 2014), fluids such as brine and partial melt (e.g., Hill et al., 2020;Wannamaker et al., 2008), and the presence of small volumetric proportions (1%-5%) of interconnected phases such as graphite or sulfides at the grain boundary scale (e.g., Hill et al., 2021;Mareschal et al., 1992;Roots et al., 2022).Thus, MT data can provide constraints on the architecture and composition of the lithosphere as well as on its thermal and geochemical state, which in turn provide constraints on its tectono-magmatic history.Here we present the results of inversion of MT data across the western Superior Craton (Figure 1), and from the generated 3D resistivity model we investigate syn-and post-orogenic processes which led to the construction and stabilization of the region.

Geological Setting
The Superior Craton of Canada forms the core of the North American Craton.The western portion of the Superior Craton (hereafter referred to as the western Superior), is comprised of largely ∼E-trending metasedimentary, gneissic-plutonic, and volcano-plutonic terranes (Figure 1), with a generally southward younging direction (Percival et al., 2006(Percival et al., , 2012)).The geodynamic processes active during the formation of the Superior Craton remain controversial, with the early stages of crustal growth being attributed to either modern-style subduction, mantleplume driven rifting, or a combination thereof (e.g., Mole et al., 2021;Percival et al., 2006;Strong et al., 2023).Nonetheless, many studies in the western Superior have resulted in a well-defined framework for construction into its current configuration via a series of late accretionary events during the terminal stages of construction of the region, each mostly followed by processes of collision, synorogenic sedimentation and magmatism, and concomitant deformation (see Percival et al., 2006Percival et al., , 2012 and references therein for further details).At ∼2,720 Ma, the northern-most subprovinces (North Caribou, Island Lake, Oxford-Stull) combined to form a composite terrane thought to be the core of a Mesoarchean proto-craton against which the southern terranes later accreted.The composite Wabigoon-Winnipeg River terrane amalgamated with this northern proto-craton during the ∼2,700 Ma Uchian orogeny, which also resulted in the trapping and deformation of the synorogenic English River sedimentary dominated subprovince.Similarly, during the ∼2,690 Ma Shebandowanian orogeny, the volcanic dominated Wawa subprovince collided with the aforementioned subprovinces to the north, while sedimentary rocks were deposited in an intervening accretionary wedge that comprises the Quetico subprovince.Regional amalgamations culminated with the collision of the Minnesota River Valley subprovince from the south, termed the Minnesotan orogeny, which occurred at ∼2,680 Ma and is the last major lithotectonic event described during the construction of the western Superior Craton.
Late to post-orogenic thermal and magmatic events resulting in emplacement of intrusive suites of both crustal and mantle derived material have been recorded throughout the western Superior (e.g., Corfu et al., 1995;Tóth et al., 2022;Whalen et al., 2004).These magmatic pulses often occur near the suture zones stitching together the various composite terranes of the western Superior.The mechanism of the generation of these melts and magmas are thought to be related to either the injection of lower-crustal and mantle derived melts due to ongoing subduction (Krogh, 1993), lithospheric inversion (Percival & Pysklywec, 2007), and/or crustal delamination and extensional collapse of locally overthickened crust (Moser et al., 1996).
Cessation of major magmatic events within the western Superior Craton occurred at ∼2,600 Ma and the craton has since remained relatively stable, with a few Proterozoic episodes of continental growth and deformation occurring Tectonics 10.1029/2023TC008110 along its margins.Deformation and further accretion along the western and southern margins of the Superior Craton occurred as result of the ∼1,800 Ma Trans-Hudson (Hoffman, 1988) and Penokean (Schulz & Cannon, 2007) orogenies, respectively.The intracratonic basins along these margins have been largely infilled by Phanerozoic sediment (Figure 1).Finally, the ∼1,100 Ma Keweenawan rifting event formed the Mid-Continent Rift, a failed intracontinental rift which resulted in significant magmatic influx whose northern arm is characterized on the surface by the Nipigon Embayment (Bedrosian, 2016).

MT Methods
The MT method is a passive source electromagnetic (EM) method that exploits the natural variations in the Earth's EM field caused by interactions between the solar wind and the Earth's magnetosphere and by global lightning activity.The resulting broadband EM waves propagate within the waveguide created by the Earth's surface and ionosphere, with partial transmission into the subsurface.The transmitted magnetic and electric fields within the Earth produce secondary fields, the amplitude and orientation of which depends on the subsurface resistivity structure.Measurement of the horizontal electric and magnetic fields, E and H respectively, allows estimation of the MT impedance tensor Z, using the relation E = ZH (Chave & Jones, 2012).The elements of Z are complex, period dependent values describing the relationship between E and H, and depend on the resistivity ρ (or conductivity, σ = 1/ρ) of the Earth.There is also a relationship between the vertical and horizontal components of the magnetic field H z = K • H, where H z is the vertical component of the magnetic field and K is the vertical magnetic transfer function (VTF), which is sensitive to lateral variations in resistivity.The phase tensor method (Caldwell et al., 2004) is commonly used to determine first order characteristics of MT data.The phase tensor, Φ, may be derived from the impedance tensor Z and displayed as ellipses.The relative size of the major and minor axes (i.e., the ellipticity) of the ellipse may indicate a preferred direction of current flow at a depth corresponding to the period of the data (Bibby et al., 2005).Secondary parameters such as Φ 2 , given by the geometric mean of the minimum and maximum phases of Φ, may be used to color the ellipses to provide an estimate of increasing (Φ 2 < 45°) or decreasing (Φ 2 > 45°) resistivity with depth.

MT Survey
The western Superior MT array used in this study consists mainly of long-period and broadband measurements made from 1997 to 2000 during the Lithoprobe project (Ferguson et al., 2005;Natural Resources Canada, 2023).The 2018 Metal Earth project collected broadband and AMT stations across the Superior Craton (Hill et al., 2021;Roots et al., 2022), including within the western Superior (Metal Earth, 2024).A small subset of the broadband stations was used to infill gaps in the Lithoprobe data coverage, resulting in a total of 185 stations across the western Superior (Figure 1).The reliability of features imaged at depth by MT inversion relies partly on the survey aperture (the lateral extents of the station coverage; Chave & Jones, 2012), and so a further 191 stations were included to the south from the Earthscope MT array (USArray, Kelbert et al., 2011Kelbert et al., , 2019) ) to ensure adequate data constraints on the deep resistivity structure, giving a total of 376 stations in the area considered.All stations include the three magnetic and two electric components, and were processed using standard robust processing techniques (Jones et al., 1989) resulting in the full impedance and VTFs.
The regions covered by the USArray MT survey are not within the scope of this project, and have been well described elsewhere (e.g., Bedrosian & Feucht, 2014;Kelbert et al., 2019;Yang et al., 2015).Furthermore, broadband MT stations in the far north of the survey area were included largely to extend the survey aperture and thus provide greater data constraints at depth within the primary areas of interest, and the regions covered by these stations are likely not spatially well resolved by the available coverage.Therefore, the discussion of the results focuses mainly on the primary area of interest (area shown in Figure 1b).

MT Inversion
Inversion of the MT data was performed using the ModEM inversion algorithm (Egbert & Kelbert, 2012;Kelbert et al., 2014).A number of inversions were performed with varying combinations of inversion parameters (e.g., model smoothing parameters), station subsets, data types (i.e., Z-only, K-only, Z + K), nominal cell size, period range, and starting model.Most of the inversions achieved reasonable levels of normalized root-mean square misfit (nRMS < 2), however, many were rejected as the modeled resistivity structures were not realistic (e.g., scale of structures was unreasonable when compared to the expected geology at various depths) and/or did not sufficiently reproduce key qualitative characteristics of the MT data (e.g., induction arrow/phase tensor orientations).Among the models that were produced, a subset were selected based on overall nRMS, the distribution of misfit among stations, periods, and data components, as well as the qualitative fit to phase tensor and induction arrow patterns for further comparison and analysis.Preference was given to those models which included data well outside the primary area of interest (i.e., the northern-most Lithoprobe stations and USArray stations to the south of the craton) in order to ensure resolution of and confidence in the large-scale structure and extents of the imaged features.This criterion also allowed comparison of the modeling results with those from previously published work (e.g., Bedrosian, 2016;Yang et al., 2015).The suite of models were compared against each other and significant similarities in the imaged resistivity structures were found in the majority of the models, suggesting that these features were robust and insensitive to changes in inversion parameters, and therefore more likely to be required by the data.A final, preferred model was selected based on overall model fit and congruence with relevant geological information, for example, alignment of the imaged resistivity structures with prominent crustal geological features.
The preferred model presented herein was generated using the nested modeling approach of ModEM, whereby an inversion mesh only slightly larger than the area covered by the MT stations is nested inside a coarser forward modeling mesh extending several skin depths away from the survey, ensuring that the MT boundary conditions are satisfied and that the calculated fields are accurate.The larger model consisted of a 154 × 119 × 73 mesh in the north-south, east-west, and vertical directions, respectively, comprised of rectilinear cells with nominal side lengths of 20 km, extending over 5,000 km in each horizontal direction.The nested inversion mesh (covering the area shown in Figure 1a) is 144 × 141 × 72 with a nominal cell size of 10 km within the primary area of interest, with cell sizes gradually increasing to 20 km within the regions to the north and south.The starting (and reference) model consisted of a 500 Ωm half-space (the average apparent resistivity of the inverted data) to 410 km depth at the mantle phase transition (Yoshino & Katsura, 2013), below which the resistivity is set to 20 Ωm.Topography was not incorporated as the changes in elevation are negligible over the large area considered, however, oceans and large bodies of seawater were incorporated in both the larger and nested meshes by setting the corresponding cells to a resistivity of 3.2 Ωm (bathymetry was downloaded from https://coastwatch.pfeg.noaa.gov).

10.1029/2023TC008110
The inversion used to generate the preferred model used 23 periods logarithmically spaced from 8 to 15,385 s for Z; for K, periods from 8 to 5,128 s were included to limit finite source field effects, which are more prevalent for K than for Z data (e.g., Yang et al., 2015).Error floors were applied using the following: for all components of Z at each frequency, ω, while a flat error floor of 0.05 was applied for K.The K data at the broadband stations is generally poor at the periods considered, and so was only included from the long-period sites.

Data Fit
The inversion was run using a sequential approach.First, a Z-only impedance was run, which converged in 131 iterations to an nRMS of 1.31.The model from the 70th iteration (nRMS of 1.7), chosen such that the modeled resistivities generally match the impedance data without overfitting, was then used to invert both Z and K data, converging in 87 iterations to an nRMS of 1.64.
The distribution of RMS misfit is spread fairly evenly among the stations (Figure 2a).The nRMS of the impedance data is 1.2 and is generally well fit across the survey.A few stations have higher than average impedance nRMS values (>3), which in general corresponds to poorly fit impedance amplitudes due to static shift effects which the coarse mesh of the model was unable to reproduce.The impedance misfit is notably higher for short periods (<20 s; Figure 2b) due to the inability to adequately recover shallow inter-station structures as a result of the mesh size and smoothness requirements of the inversion.Nonetheless, phase tensor ellipses for the observed data are reproduced by the model (Figure 3), suggesting that the primary characteristics of the data are well matched.
The misfit of the VTF data (Figure 2, green line) is typically higher than that of the impedance data, with an nRMS of 1.98.Whereas the impedance fit is distributed relatively equally among all stations, the VTF misfit is higher for those stations to the north of the area of interest (Figure 1b), likely due to the combination of smaller station spacings and increased complexity of the structures relative to those within the western, central, and eastern blocks.Nonetheless, the induction arrow directions are generally well matched (Figure 3; see also Figure S3 in Supporting Information S1), indicating that the geometries of the large-scale resistivity structures fit the data.As with the impedances, the nRMS of the VTF data are higher at the short periods.Additionally, there is a minor increase in misfit at the longest periods where induction arrows largely point to the west (Figure 3b), and therefore the increase in misfit is attributed to conductive structures which lie outside the modeled area.

Seismic Reflection Data
Lithoprobe seismic reflection data were collected in 1996-1997 along the ∼600 km long "Line 1" (Figure 1b, red dashed line).The data were collected to a two-way travel-time of 16 s and processed using standard techniques (Natural Resources Canada, 2019;van der Velden & Cook, 2005;White et al., 2003) resulting in reflection images extending into the uppermost mantle (∼50 km depth).Recent re-processing of the data along Line 1 was performed to remove artifacts from the crooked acquisition geometries to facilitate comparisons of the reflection images with those from modern seismic attribute analysis and subsequent interpretation of the polyphase deformation history of the western Superior (Calvert et al., 2023).

Inverse Model
A section showing the depth-integrated resistivity of the modeled area (excluding the larger model within which the nested inversion was performed) in the upper-to mid-crust (<20 km depth) is shown in Figure 4, while Figure 5 presents depth slices from the mid-crust to the lithospheric mantle, restricted to the primary region of interest (Figure 1b).The following sections discuss the modeling results at different depth intervals, with the exception of the C1 feature which is discussed first.

Upper Crust
Due to the relatively large station spacing, cell sizes, and long periods used in the inversion, the details of the upper crust resistivity structure at shallow depths (<20 km) are not resolvable by the available data.However, the calculated thickness-weighted integrated resistivity section (i.e., the inverse of the conductance divided by the thickness of the integrated section; Figure 4) reveals first-order characteristics consistent with those imaged or inferred by previous works (e.g., Bedrosian, 2016;Ferguson et al., 2005;Yang et al., 2015), which lends some confidence that the results presented here are reliable.
The craton boundaries are well resolved (up to the limitations of the station spacing) with the upper crust of the western Superior Craton having high resistivities throughout (>1,000 Ωm; Figure 4).Small zones of lower resistivity on the order a few model cells are seen throughout the model.Many of these low resistivity zones are spatially correlated to the greenstone belts of the Wabigoon subprovince in the central block and are likely a result of galvanic distortion resulting from small-scale heterogeneities within these geologically complex regions that are below the resolving power of the available data.Along the western edge of the craton, resistivities drop by an

Lower Crust
Coherent large-scale resistivity patterns begin appearing in the model in the mid-to lower crust (20-45 km; Figures 5a-5c).The lower crust of the central block is characterized by curvilinear high and low resistivity zones, generally sub-parallel to the subprovince boundaries.In particular, two significant >250 km long low resistivity (1-100 Ωm) anomalies, C2 and C3, with intervening high resistivity zones, R2 and R3, are imaged (Figures 5b  and 5c).C2 forms a U-shaped low resistivity feature with a resistive core (R1) that underlies the English River and Winnipeg River subprovinces.The northern arm of C2 extends for ∼200 km along the English River subprovince, before bending south and then extending west along the boundary between the Winnipeg River and Wabigoon subprovinces.The vertex of C2 roughly coincides with where the English River and Winnipeg River subprovinces pinch northwards toward the Uchi subprovince.The C3 low resistivity feature is sub-parallel and spatially correlated with the northern edge of the Quetico subprovince and is flanked on the north and south by high resistivity (>10,000 Ωm) features, R2 and R3, respectively.R2 runs sub-parallel to C3 until ∼100 km west of the Nipigon Embayment, where it terminates south of the vertex of C2, whereas R3 continues along the entire length of C2, underlying most of the Quetico subprovince within the western Superior.
The mid-to lower crust of the northern and southern blocks are relatively homogenous, having generally high resistivities (>10,000 Ωm) with only a few isolated regions of lower resistivity (Figures 5a-5c).The lower crust of the eastern block is more heterogenous, with reduced resistivities beneath the Nipigon Embayment corresponding to C3 and a set of low resistivity zones (C4) imaged at the eastern extent of the model.While indications of C4 are visible in the induction arrow data, it is largely unconstrained by the available stations and so will not be discussed further.

Lithospheric Mantle
The orientations of the imaged resistivity structures at mantle depths are markedly different from those within the crust.Where the resistivity structures within the crust are approximately east-trending and sub-parallel to the subprovince boundaries, those within the mantle are oriented roughly north-south, with the upper mantle model depths (63-89 km; Figures 5d and 5e) acting as a transition zone, likely due to the regularization of the inversion algorithm.At 89 km depth (Figure 5d), the features imaged in the lower crust transition into a series of approximately north-south alternating high and low resistivity bands (C5a,b,c).The low resistivity bands are 60-100 km wide with resistivities of 1-100 Ωm and are largely contained within the central block with each terminating beneath the Quetico subprovince to the south and the Uchi subprovince to the north, with the exception of C5a which extends into the North Caribou subprovince.
With the exception of regions where the C5 conductors and associated resistive bands extend, the resistivities imaged at mantle depths in the model are relatively simple compared to those in the central block.In the western block, the C1 conductor extends along the western edge of the craton to a depth of 129 km.Beneath the MT stations in the western block the model resistivities gradually decrease from >10,000 Ωm at 63 km to <1,000 Ωm at 244 km.Similarly, with the exception of the northern-most extents of C5a, the majority of the northern block exhibits high resistivities (>10,000 Ωm) at depths of 63-89 km and notably lower resistivities (100-500 Ωm) at 244 km.The southern block has similar characteristics, with high resistivities (5,000-10,000 Ωm) at 63 km, decreasing to 30-300 Ωm at 244 km, with lower resistivities beneath the Wabigoon and Quetico subprovinces.In the eastern block, the C5c conductive band beneath the Nipigon Embayment has low resistivities (<500 Ωm) throughout the mantle.The low resistivities of C4 extend from 63 to 89 km depth; beyond this depth, the region east of the Nipigon Embayment exhibits moderate to high resistivities (1,000-10,000 Ωm).
The imaged resistivity structure of the lithospheric mantle beneath the western Superior is highly complex.The alternating high/low resistivity banding (Figures 5f and 5g) has few plausible geological explanations, and is likely an artifact of electrical anisotropy within isotropic models (Wannamaker, 2005).Additional work is required on the investigation of the MT data at these depths and will be the topic of a future study.

Seismic Reflection Data
Lithoprobe seismic Line 1 extends ∼600 km between the Wawa and North Caribou subprovinces, intersecting many of the resistivity features described above.The migrated section and interpretive lines (Calvert et al., 2023) are presented and overlain on the coincident depth slice through the MT resistivity model (Figure 6).The seismic

Discussion
The resistivity model reveals several notable characteristics of the crust of the western Superior.A clear distinction exists between the subprovinces forming the cratonic core (i.e., the northern block) and those forming the belt-like subprovinces to the south.Though the finer details of the resistivity structure of the northern block are unresolved due to the large station spacing, the model reveals high resistivities (>10,000 Ωm) and relatively minor lateral variation throughout the crust, similar to previous interpretations of the MT data (Ferguson et al., 2005).The high resistivities observed within the crust are consistent with those expected for the composition of felsic intrusive suites interspersed with relatively small supracrustal belts (Percival et al., 2012).A few zones of lower resistivity (<1,000 Ωm) exist; however, these generally span no more than a single station, and are therefore not sufficiently constrained by the available data.
The relative homogeneity of the craton core in the north is juxtaposed against the more complex resistivity structure of the terranes to the south, which record a history of crustal growth within the western Superior.The belt-like subprovinces within the central block are characterized by bands of high (>10,000 Ωm) and low (<100 Ωm) resistivity within the mid-to lower crust (i.e., C2, C3, R1-R3; Figures 5b and 5c) which trend subparallel to and beneath or slightly offset from the subprovince boundaries.Alternating bands of high/low resistivity similar to those imaged in the lower crust (Figure 5c) may be associated with electrical anisotropy which can arise from aqueous fluids, melts, and/or solid-state conductors such as graphite emplaced within preferentially oriented structures such as dikes or shear zones (e.g., Bedrosian et al., 2019;Wannamaker, 2005).However, thermal stability of the region since the Neoarchean (Percival et al., 2012) precludes the presence of melt, and free aqueous fluids are unstable at lower-crustal depths over the relevant time scale (Rudnick & Fountain, 1995;Yardley et al., 2010).Furthermore, electrical anisotropy with a WSW-ENE orientation and large spatial extent exceeding 150,000 km 2 would require a vast network of near-vertical conductors (e.g., graphitized shear zones) with WSW-ENE orientations, which is inconsistent with the structural fabric of the region (Percival et al., 2012 and references therein) and the pervasive sub-horizontal lower-crustal reflections observed in the seismic data (Figure 6; Calvert et al., 2004Calvert et al., , 2023;;White et al., 2003;van der Velden & Cook, 2005).As such, the observed east-west trending lower-crustal resistivity structures are interpreted as reflecting real isotropic features rather than artifacts of anisotropy.
Conductive anomalies within the lower crust and upper mantle of Archean terranes are commonly attributed to domains containing interconnected conductive mineral phases (e.g., graphite, sulphides), perhaps limited to grain-boundaries along sub-horizontal shear zones (e.g., Adetunji et al., 2023;Evans et al., 2011;Haugaard et al., 2021;Heinson et al., 2021;Hill et al., 2021;Roots et al., 2022;Selway, 2014), and are the most likely cause of the C2 and C3 low resistivity features.The mechanism by which the conductive phases were emplaced is more enigmatic, though additional geological and geophysical data can be used to constrain potential mechanisms.
Lithoprobe seismic reflection data imaged north-dipping reflections at Moho depths beneath the western Superior Craton that have been interpreted to represent either accretionary suture zones and/or syn-to post-accretionary extensional shear zones at the lower crust to upper mantle boundary (Calvert et al., 2004;White et al., 2003).
Ductile flow and extension along associated mid-to lower-crustal shear zones is inferred from sub-horizontal to listric reflections, indicative of pervasive strain at these crustal levels (Calvert et al., 2023), and which overlap with the C2 and C3 low resistivity features (Figure 6).Modeling of gravity data across the Uchi and English River subprovinces coincident with the western portion of the C2 low resistivity feature revealed undulations in Moho depth as well as a distinct lack of intracrustal mass anomalies providing evidence for a syn-or post-accretionary thermal weakening of the crust and subsequent lateral ductile flow (Calvert et al., 2004;Nitescu et al., 2006).This thermal weakening may be linked to increased heat flow and the emplacement of igneous rocks generated through partial melting of the lower crust triggered by crustal delamination (Corfu et al., 1995;Stevenson et al., 1999).Late-to post-orogenic delamination of locally overthickened crust provides a viable mechanism for thermal weakening, extension, and differentiation of the crust (Meissner & Mooney, 1998;Rey, 1993) and while delamination is less common in modern orogens, it may have been the predominant form of lithospheric recycling in the Archean (Chowdhury et al., 2017;Johnson et al., 2014).Delamination of crustal material results in upwelling mantle and melting of the lower crust (e.g., Chowdhury et al., 2017), which would preferentially partition volatile components (e.g., CO 2 and H 2 S) into the melt/fluid phase (Hill et al., 2021;Krogh, 1993).These carbon and sulfur rich fluids could then migrate upwards through the crustal column, likely through conduits defined by suture zones and associated crustal shear zones (e.g., Figure 7a), precipitating graphitic carbon and possibly sulphides upon passage through low fO 2 crustal rocks (Fu & Touret, 2014;Luque et al., 2014).Reactivation and extension along the lower-crustal shear zones as a result of thermal weakening may further enhance interconnection of conductive phases (Figure 7b).While not evident in the presented model, geological studies (e.g., Tóth et al., 2022) suggest related melts may be linked to late stitching plutons along terrane boundaries.Subduction and slab breakoff may provide an alternative geodynamic model consistent with the observed geophysical signatures (e.g., McGary et al., 2014); however, geodynamic modeling studies suggest shallow (40-80 km depth) slab necking and breakoff is required for the generation of mantle melts, which may be uncommon except for very young (<30 Ma) plates (Duretz et al., 2011;Freeburn et al., 2017), and therefore delamination is the preferred interpretation.
The available geological and geophysical data suggest a similar history to the south within the Wawa-Abitibi, Quetico, and Wabigoon subprovinces, coincident with the R2, C3, and R3 resistivity features (Figures 5b, 5c, and  6).The low resistivity C3 anomaly extends from mid-crust to upper mantle depths, running subparallel to the northern border of the Quetico subprovince for ∼600 km.The continuity of C3 suggests that it is largely the result of Neoarchean history of the Quetico and its proximal domains, and that later events (e.g., the Mid-Continent Rift) did not significantly modify the corresponding large-scale geoelectric signatures.It is possible that the enhanced conductivity of C3 may be related to the collision between the Quetico and domains to the north.Burial and metamorphism of carbon-rich sedimentary rocks to mid-crustal depths would result in abundant graphite-bearing rocks (e.g., Cherevatova et al., 2014;Heinson et al., 2022); however, sub-horizontal reflections beneath the basin (Figure 6) suggests that the Quetico basin is relatively shallow (∼12 km deep; White et al., 2003) compared to the depth of C3, which would require significant tectonic displacement by underthrusting of Quetico sedimentary rocks beneath the Wabigoon and Marmion subprovinces.Therefore, while underthrust sediments may provide a minor source of carbon in the lower crust, it is not likely the sole cause of the C3 anomaly.Intrusions within the Quetico subprovince as well as within the neighboring Wabigoon and Abitibi-Wawa subprovinces display characteristics of mantle derived magmas (e.g., Stevenson et al., 1999;Tóth et al., 2022;Whalen et al., 2004).Specifically, ∼2,690 Ma sanukitoid suites along the border between the Quetico and Wawa-Abitibi subprovinces suggest magmatism related to mantle upwelling along and north of the inferred suture beneath the Quetico (Percival et al., 2006;Tóth et al., 2022;Whalen et al., 2004).Seismic reflection data beneath the Quetico and Wabigoon subprovinces image pervasive shallowly dipping reflections (Figure 6; see also White et al., 2003;van der Velden & Cook, 2005), suggestive of mid-to lower-crustal shear zones activated (or reactivated) during extension in response to late to post-orogenic magmatism (e.g., Calvert & Doublier, 2018;Hill et al., 2021).Such extension/shearing in addition to an influx of mantle-derived CO 2 -rich fluids and possibly release of additional CO 2 via retrograde metamorphism of underthrust sediments (e.g., Heinson et al., 2021) provides a mechanism for the generation and widespread interconnection of graphite-bearing films in the mid-to lower crust within an otherwise resistive Archean crust.As with the C2 feature to the north, we suggest that the generation and emplacement of mantle-derived melts and subsequent crustal extension may have been triggered by crustal delamination or shallow slab breakoff.
The crustal and upper mantle resistivity structure of the western Superior Craton is indicative of the late stages of Archean crustal growth and cratonization via accretionary processes with late episodic magmatism related to partial melting of the lower crust and/or direct input of mantle derived melt (Figure 7a).Subsequent extension and reactivation of lower-crustal shear zones resulting in ductile flow of the mid-to lower crust may have served to increase interconnection of conducting phases (Figure 7b), resulting in the linear low resistivity features imaged herein (C2, C3; Figures 5b and 5c).We speculate that these syn-to post-orogenic events may have been triggered by collapse and delamination of locally overthickened crust and/or gravitational instability within a hotter Archean geothermal regime (e.g., Hill et al., 2021;Johnson et al., 2014).Furthermore, the imaged resistivity features span a large proportion of the western Superior Craton, and thus we suggest that the relevant geodynamic mechanisms represent craton-scale phenomenon are fundamental to models of Archean crustal growth and stabilization.While the model and data set presented herein lack sufficient resolution in the upper crust, the occurrence of abundant late hydrous granitoids and sanukatoids with mantle affinity in the study area (e.g., Tóth et al., 2022) suggests that the proposed mechanism could also result in progressive migration of enriched fluids into the upper crust.However, the apparent correlation between significant deposits and crustal conductors within the western Superior is noticeably weaker than that found elsewhere (Jørgensen et al., 2022;Kirkby et al., 2022).Thus, questions remain as to whether the apparent lack of metal endowment within the western Superior Craton is a result of insufficiently fertile melts, a lack of focusing mechanisms (e.g., melts spread across a large area rather than focused into economically viable mineral deposits), or due to relatively limited exploration of the region.Identification of corresponding upper-crustal geophysical signatures (or lack thereof) within the western Superior, and in particular within the footprints of the C2 and C3 anomalies, could provide a link between the proposed geodynamic model and mineralization potential of the western Superior, and would therefore have significant implications on the understanding of the mineral prospectivity of the region.

Figure 1 .
Figure 1.Simplified geological map of the study area including (a) all inverted MT stations and (b) the primary area of interest as discussed in the text.Solid white line indicates Canada-US border.White dashed lines and corresponding white text denotes blocks of stations discussed in the text: CB-Central Block; EB-Eastern Block; NB-Northern Block; SB-Southern Block; WB-Western Block.Western Superior subprovince boundaries are shown by solid black lines: ER-English River; IL-Island Lake; Mm-Marmion; NC-North Caribou; OS-Oxford-Stull; Qt-Quetico; Uc-Uchi; WBG-Wabigoon; Ww-Wawa; WR-Winnipeg River.Red text corresponds to regions affected by Proterozoic tectono-magmatic events; MCR-Mid-continent Rift; NE-Nipigon Embayment; THO-Trans-Hudson Orogen.Exterior coordinates projected to Canada Atlas Lambert (EPSG: 3979).Map data modified from the compilation of Montsion et al. (2018).
striking linear feature C1 is imaged along the western edge of the model, with resistivities <1 Ωm and extending from depths between 31 and 129 km.This feature is consistent with the large amplitude, westward oriented induction arrows seen at stations within the western block (Figure3), however, the depth extents of C1 raise questions as to its validity.Previous MT surveys and 2D inversions of data within the Trans-Hudson Orogen have imaged extremely low resistivity anomalies (<1 Ωm) within the crust and mantle associated with the North American Central Plains anomaly (e.g.,Jones et al., 2005), as well as within the mid-to lower-crust of the Thompson belt (e.g.,Gowan et al., 2009).Gowan et al. (2009) also imaged a low to moderate resistivity (10-100 Ωm) feature in the crust/upper mantle along the western edge of the Superior craton, however, no feature with the extremely low resistivities and depth extents of C1 has previously been imaged.While it is possible that C1 represents a previously undiscovered low resistivity anomaly, it is more likely that C1 represents the projection of multiple more distal low resistivity structures into the considered model space.

Figure 3 .
Figure 3. Phase tensor and induction arrow plots at periods of 40 and 420 s for the observed data (a, b) and modeled responses (c, d).Ellipses are colored by Φ 2 .

Figure 2 .
Figure 2. Normalized root mean square misfit statistics for the final inverse model organized by (a) geographic location of the stations from south to north and (b) by period.

Figure 4 .
Figure 4. Thickness-weighted integrated resistivity of the top 20 km of the preferred inverse model.The area presented corresponds to that in Figure 1a; the white dashed box corresponds to the area in Figure 1b.Solid black lines correspond to major terrane boundaries including subprovince boundaries, the craton edge, and the surface expression of the Mid-Continent Rift.

Figure 5 .
Figure 5. Plan view slices at representative depths through the preferred inverse model at depths of (a) 20 km, (b) 31 km, (c) 45 km, (d) 63 km, (e) 89 km, (f) 129 km, (g) 189 km, and (h) 244 km.White dashed lines in (a) correspond to the blocks of stations as discussed in the text and in Figure 1b; red dashed line corresponds to Lithoprobe seismic line 1.White dashed lines in (b)-(h) highlight the features discussed in the text.Text in (h) corresponds to subprovince boundaries as in Figure 1.Color scale is the same as that in Figure3.The area presented in these slices corresponds to that in Figure1b.
described and interpreted in detail elsewhere (e.g.,van der Velden & Cook, 2005;White et al., 2003) and so we limit our discussion to features which are relevant to the interpretation of the resistivity volume.The profile exhibits relatively low reflectivity within the upper crust (0-10 km); where reflections are present, they are generally short and discontinuous.Reflectivity increases in the mid-crust (10-25 km) with dense packages of largely sub-horizontal to shallowly north dipping reflections evident throughout much of the section.Reflections decrease slightly in dip with increasing depth.Reflections beneath the Quetico, Marmion, and Winnipeg River subprovinces gradually become listric into the lower crust and flatten out near the Moho.A dense package of listric reflections overlap with the C3 conductive anomaly in the lower crust beneath the Marmion subprovince, while a similar set of reflections overlaps the southern portion of C2 beneath the Wabigoon/ Winnipeg River subprovinces.Reflections flatten toward the Moho, which is characterized by a variably sharp to moderately diffuse transition from high to low reflectivity at depths of 38-45 km.Beneath the Wabigoon/ Winnipeg River and Uchi subprovinces, north dipping reflections extend beneath the Moho and into the upper mantle.

Figure 6 .
Figure 6.South to north vertical slice through the 3D resistivity volume along and overlain by the Lithoprobe Line 1 seismic reflection data (no vertical exaggeration).Simplified surficial geology coincident with the section are included along the top with colors representing lithologies as in Figure 1.Overlain interpretive lines are modified from Calvert et al. (2023).

Figure 7 .
Figure 7. Schematic diagram illustrating how syn-to post-orogenic processes may have produced the observed crustal resistivity signature of the western Superior Craton.(a) Shallow slab breakoff or delamination of overthickened crust results in melting of the lower crust and/or upper mantle.Incompatible elements (including CO 2 and possibly H 2 S) are preferentially incorporated into the resulting melt phase, and migrate upwards along zones of weakness, likely corresponding to crustal shear zones formed during amalgamation of the subprovinces.(b) Increased heat flow at the crust-mantle boundary results in thermal weakening of the mid-to lower crust and extension along reactivated crustal shear zones, promoting formation and interconnection of grain boundary conductive films.Upwards migration of melts results in late stitching plutons along terrane boundaries.
(Kirkby et al., 2019;Krieger & Peacock, 2014)iga Geoscience, under contract to Complete MT Solutions.Much of the Lithoprobe MT data collection was completed by Phoenix Geophysics.We would like to thank and acknowledge the Digital Research Alliance of Canada for high-performance computational support.We acknowledge the use of MTPy(Kirkby et al., 2019;Krieger & Peacock, 2014)for conversion of the data from spectral to impedance estimates.We would also like to thank and acknowledge C. Farquharson and an anonymous reviewer for their helpful comments and reviews, as well as M. Rushmore for editorial handling of this manuscript.The Metal Earth project is supported through the Canadian government via the Canada First Research Excellence Fund.This is Metal Earth