Genesis of Sublayer, Footwall Breccia, and Associated Ni-Cu-Platinum Group Element Mineralization in the Sudbury Igneous Complex

The Sudbury structure is one of the oldest (~1.85 Ga), largest (up to 260-km diam), and best-exposed terrestrial impact structures (Dietz, 1964; Grieve and Therriault, 2000; Spray et al., 2004) and contains some of the world’s largest magmatic Ni‐Cu‐platinum group element (PGE) sulfide deposits (e.g., Naldrett, 2004; Lightfoot, 2016). It comprises the Sudbury Igneous Complex, overlying fallback breccias, suevites, and basin‐fill sediments of the Whitewater Group (e.g., Muir and Peredery, 1984; Ames et al., 2002), and underlying pseudotachylitic Sudbury breccias, anatectic footwall breccia, and contact‐metamorphosed country‐rock breccias (e.g., Dressler, 1984; McCormick et al., 2002). The Sudbury Igneous Complex straddles the boundary between the 2.7–2.6 Ga Superior province to the north, east, and west (North Range) and the 2.5 to 2.2 Ga Southern province to the south (South Range) (Fig. 1A). The Superior province in this area consists mainly of felsic‐mafic gneisses of the 2.7 to 2.6 Ga Levack Gneiss Complex (e.g., Krogh et al., 1984), the 2.6 Ga Cartier granite (Meldrum et al., 1997), and minor mafic‐ultramafic rafts/intrusions (e.g., Moore et al., 1993, 1994, 1995). The Southern province comprises metabasalts, metarhyolites, and metasedimentary rocks of the 2.5 to 2.2 Ga Huronian Supergroup (e.g., Young et al., 2001), 2.5 Ga felsic granitoids (e.g., Dressler, 1982), and mafic‐ultramafic intrusions of the 2.5 Ga East Bull Lake (e.g., James et al., 2002) and 2.2 Ga Nipissing (e.g., Lightfoot and Naldrett, 1996) Suites.

The Sudbury Igneous Complex comprises a main mass, the underlying discontinuous inclusion‐rich sublayer and footwall breccia, and radial and concentric offset dikes (Fig. 1A). The main mass is differentiated into a lower noritic unit, a transitional quartz gabbroic unit, and an upper granophyric unit and has been interpreted to be the remnant of a crystallized impact melt sheet (e.g., Dietz, 1964; Grieve et al., 1977). The sublayer is an up to 700-m-thick, heterogeneous, inclusion‐rich, fine‐ to medium‐grained noritic to gabbronoritic unit that is localized in embayments, funnels, and troughs along the basal contact of the main mass (see review by Lightfoot, 2016) (Fig. 1B). Offset dikes contain an analogous inclusion-rich quartz diorite (simplified hereafter to “inclusion quartz diorite”) core with inclusion‐poor quartz diorite (simplified hereafter to “quartz diorite”) margins. The sublayer on the North Range and parts of the South Range is underlain by locally mobilized anatectic footwall breccia (e.g., Coats and Snajdr, 1984; McCormick et al., 2002). The sublayer, footwall breccia, and inclusion quartz diorite all contain variable amounts of Fe‐Ni‐Cu sulfides. On the North Range, the sublayer almost always overlies footwall breccia and associated clastic and semimassive contact ores, whereas on the South Range, the sublayer sometimes overlies the footwall breccia but often grades downward directly into inclusion-rich semimassive to massive contact ores (e.g., Souch et al., 1969; Naldrett et al., 1984). Offset ores range from coarse blebby disseminated sulfides to inclusion‐rich semimassive sulfides (see review by Lightfoot, 2016).

The abundant inclusions within the sublayer, footwall breccia, and inclusion quartz diorite comprise gneissic, diabasic, metavolcanic, and metasedimentary xenoliths derived from local country rocks (e.g., Souch et al., 1969; Grant and Bite, 1984) and a variety of mafic‐ultramafic inclusions. The latter comprise harzburgite, phlogopite lherzolite, feldspar lherzolite, wehrlite, olivine clinopyroxenite, orthopyroxenite, olivine melanorite, olivine melagabbronorite, and olivine gabbro on the North Range (Rae, 1975; Wang, 2019; Wang et al., 2020) and harzburgite, olivine orthopyroxenite, orthopyroxenite, troctolite, and olivine melanorite on the South Range (Scribbins et al., 1984). The mafic‐ultramafic inclusions vary from a few centimeters to tens of meters in diameter (Naldrett et al., 1984; Wang, 2019) and are interpreted to have multiple origins, including (1) anteliths derived from a disrupted early border phase crystallized from a mixture of impact melt and a more mafic melt, probably derived by melting of ultramafic footwall rocks (e.g., Prevec, 2000; Prevec et al., 2000,), (2) local xenoliths incorporated from layered mafic‐ultramafic intrusions in the country rocks (e.g., footwall of Levack and Bowell embayments: Card and Pattison, 1973; Pattison, 1979; Moore et al., 1993, 1994, 1995; Wang et al., 2020), and (3) exotic xenoliths from unexposed mafic‐ultramafic rocks at upper‐middle crustal levels (Rae, 1975; Scribbins, 1978; Wang et al., 2018, 2020).

There is no consensus on the mechanism(s) of formation of the sublayer. Some workers have suggested that it represents a mixture of splash‐emplaced impact melt and fragments driven up the walls of the crater (Pattison, 1979; Morris, 1984), some have argued that it crystallized from a hybrid melt formed by mixing between the impact melt and a mantle‐derived high‐ Mg magma (e.g., Naldrett and Kullerud, 1967; Naldrett et al., 1984; Zhou et al., 1997), and some have argued that it crystallized from a mixture of basal impact melt and a basaltic melt derived by melting mafic footwall rocks (e.g., Golightly, 1994; Lightfoot et al., 1997; Prevec, 2000; Prevec et al., 2000; Prevec and Cawthorn, 2002).

This study builds on work by Wang et al. (2018) on shock metamorphic features in the mafic‐ultramafic inclusions and by Wang et al. (2020) on the mineralogy and geochemistry of the mafic‐ultramafic inclusions, and focuses on the relationship between sublayer, footwall breccia, inclusion quartz diorite, and associated Ni‐Cu‐PGE sulfide mineralization within the context of the overall evolution of the Sudbury structure.

The principal characteristics of the sublayer, footwall breccia, and inclusion quartz diorite are summarized in Table 1. The key points are as follows:

1. The sublayer, inclusion quartz diorite, and footwall breccia may all be barren but are normally at least weakly mineralized.

2. All contain a wide variety of inclusions, including abundant local felsic-intermediate-mafic xenoliths; some contain olivine melagabbro anteliths; and all contain exotic mafic-ultramafic xenoliths, especially where mineralized.

3. The sublayer contains cumulus orthopyroxene and plagioclase with ophitic to locally porphyritic textures; inclusion quartz diorite contains noncumulus orthopyroxene and plagioclase sometimes with radiating, locally spherulitic, and rarely quench textures; and footwall breccia has a felsic-intermediate anatectic matrix.

4. The footwall breccia, the Frood-Stobie concentric breccia belt, and to a lesser degree radial offset inclusion quartz diorite all contain some sulfide clasts, indicating reworking of previously crystallized sulfides.

The sublayer and footwall breccia are localized within depressions at the base of the main mass, some of which are broad shallow embayments (e.g., Levack, Trill), some of which are narrow deep troughs (e.g., Creighton), and some of which are funnels to offset dikes (e.g., Copper Cliff, Foy, Whistle, Worthington). The thicknesses of the sublayer and footwall breccia are highly variable and are related to topography along the footwall contact. They are absent along relatively planar segments of the contact, 10 to 100 m thick within shallow embayments (e.g., Little Stobie, Murray), and up to 700 m thick within deep embayments (e.g., Trill, Whistle). The geologic features of major sublayer-hosting structures on the North and South Ranges are summarized in Appendix Table A1.

The inclusion quartz diorite and mobilized footwall breccia (also known as “metabeccia”; e.g., Lafrance et al., 2014) are localized within the proximal to medial parts of offset dikes. Most radial dikes have inclusion- and sulfide-poor quartz diorite margins and inclusion ± sulfide-rich inclusion quartz diorite cores, although inclusion quartz diorite sometimes occurs along the margins of the dikes, especially near fault dislocations (e.g., Grant and Bite, 1984; Lightfoot, 2016). The contacts between quartz diorite and inclusion quartz diorite are often sharp with inclusions of quartz diorite within inclusion quartz diorite (e.g., parts of Foy, parts of Copper Cliff, Trill, Worthington) but are sometimes gradational (e.g., other parts of Copper Cliff and Foy) (e.g., Grant and Bite, 1984; Lightfoot and Farrow, 2002; Lightfoot, 2016; Pilles et al., 2018).

Petrographic and geochemical characteristics of the igneous‐textured sublayer matrix have been described by Wang et al. (2020), and an overview is provided in the Appendix, which includes Tables A1 through A4 and Figures A1 through A9. The majority of the sublayer matrix has relatively high S contents (0.2–8%) and exhibits a strong positive correlation between Ni and S, indicating that the Ni content is controlled by pentlandite. Mafic‐ultramafic inclusions contain few if any sulfides (Ni < 1,000 ppm, S < 0.5%: App. Fig. A2) (e.g., Naldrett and Kullerud, 1967; Farrell, 1997), except for melanorite and olivine melanorite anteliths in the Whistle embayment, which contain disseminated or net‐textured Fe‐rich sulfides (200–2,000 ppm Ni, 0.5–3% S: App. Fig. A2), and local gabbro xenoliths in the Little Stobie mine, which contain fine disseminations and/or fracture fillings of sulfide (Davis, 1984).

There are two end‐member models for the formation of the sulfides in sublayer, inclusion quartz diorite, and footwall breccia, parts of which are similar and parts of which are mutually exclusive:

1. Exsolution from impact melt (Fig. 2, model A): This model involves dissolution of S from target rocks during formation of superheated impact melt, exsolution of immiscible sulfide droplets during cooling (Keays and Lightfoot, 2004; Li and Ripley, 2005), settling of sulfides and inclusions to the contact, and migration along the contact into topographic embayments driven by convective currents (e.g., Lightfoot et al., 2001; Keays and Lightfoot, 2004; Zieg and Marsh, 2005) and/or gravity‐driven density flow (this study). In this model the sulfides in contact‐footwall and offset ore systems are derived from the melt sheet and are modified by incorporating some S and metals from country rocks. Upward-decreasing metal contents in the norites are consistent with removal of sulfides (Keays and Lightfoot, 2004), but the metal tenors of coarse blebby sulfides in the lower norites overlap those in the sublayer and accepted Sudbury averages (~5% Ni, ~5% Cu: Keays and Lightfoot, 2004), so the metal depletion in the norites must have occurred after formation of the majority of the ores.

2. Xenomelts from country rocks (Fig. 2, model B): This model involves volatilization of significant amounts of Pb (O’Sullivan et al., 2016; Kenny et al., 2017; McNamara et al., 2017), Sb (O’Sullivan et al., 2016), Zn‐Cd‐Rb‐Cs (Kamber and Schoenberg, 2020), and probably also much/most of the S‐Se‐Bi and significant amounts of Ag‐Cu‐Au‐ As (Lesher, 2019a) during impact, significant postimpact thermomechanical erosion of footwall rocks by super-heated impact melt (Prevec and Cawthorn, 2002; Gregory, 2005), assimilation of miscible silicate components (see Lesher and Campbell, 1993; Lesher and Burnham, 2001), and generation of immiscible sulfide xenomelts that were upgraded by reaction with the overlying melt sheet (Lesher, 2019b). In this model the sparse sulfides in the lower main mass noritic rocks and upper parts of the sublayer segregated from the overlying melt sheet, but the majority of the sulfides in the sublayer and footwall breccia were derived from eroded footwall rocks (S‐bearing Huronian‐Superior lithologies and Ni‐Cu‐PGE–rich mineralization in Nipissing and East Bull Lake Intrusive Suite intrusions).

These are end-member processes that must have occurred to some degree in any case, so we are not able to invalidate one or the other, but we will consider all of the isotopic, physical (impact), thermal, fluid dynamic, and temporal constraints on these processes and the formation of the sublayer, inclusion quartz diorite, footwall breccia, and associated Ni-Cu-PGE mineralization below.

The Sm‐Nd isotope compositions of the Sudbury Igneous Complex lithologies, mafic‐ultramafic inclusions in the sublayer and inclusion quartz diorite, mafic‐ultramafic intrusions in the footwall, and mafic footwall lithologies are plotted as ε_{Nd(1850 Ma)} values in Appendix Figure A3 and summarized in the Appendix. The data confirm a strong continental crustal signature for all components of the Sudbury Igneous Complex (e.g., Faggart et al., 1985; Dickin et al., 1992; Prevec et al., 2000), including the mafic‐ultramafic inclusions. The two Mc-Creedy West sublayer samples with near‐zero ε_{Nd(1850 Ma)} values represent a more primitive component (Prevec et al., 2000).

The Re‐Os isotope compositions of sulfides in the sublayer, mafic‐ultramafic inclusions, and some country rocks are plotted as γ_{Os(1850 Ma)} values in Appendix Figure A4 and summarized in the Appendix. The moderately to highly radiogenic γ_{Os(1850 Ma)} values of the sulfides in the sublayer indicate that Os was derived largely from old crust, not mantle or juvenile crust (Walker et al., 1991, 1994; Dickin et al., 1992). Small but unambiguous variations in Os isotope compositions of sulfides from different mines and different parts of the same mines have been attributed to assimilation of local wall rocks with different Os isotope compositions (Walker et al., 1991, 1994; Cohen et al., 2000) or to fractional crystallization of monosulfide solid solution (MSS; Morgan et al., 2002). Sulfur and Pb isotope data are not available for the same samples, but some ores have anomalous S (e.g., Crean Hill contact, Levack contact, Levack footwall, Lindsley, Victor 42N: Ripley et al., 2015) and Pb‐Pb (e.g., Blezard footwall, Morrison footwall: McNamara et al., 2017) isotope ratios, which should not change during fractional crystallization. As a result, it seems likely that most of the variations in Re‐Os isotopes also resulted from local contamination.

Returning to the question of the source, Re is less compatible than Os during mantle melting (Walker and Morgan, 1989). Subcontinental lithospheric mantle that has been depleted by substantial fractions of partial melt will consequently evolve to negative γ_Os values (Walker and Morgan, 1989; Carlson and Irving, 1994). As a result, the negative γ_{Os(1850 Ma)} values of the S‐rich gabbro and S‐poor melanorite inclusions in the Whistle mine are compatible with a parental magma derived from a Re‐poor lithospheric mantle source (see e.g., Ellam et al., 1993; Lambert et al., 1994) or from older crustal rocks with low Re/Os ratios and unradiogenic Os isotope signatures (see Cohen et al., 2000). However, the highly unradiogenic Nd signatures are not compatible with a direct contribution from subcontinental lithospheric mantle. Instead, older crustal rocks, e.g., Archean volcanic and subvolcanic rocks in this region, which have been proposed to be derived from subduction‐metasomatized mantle (Jolly, 1987; Jolly et al., 1992; Ketchum et al., 2013), must have been sampled during the impact.

The Pb-Pb isotope compositions of main mass lithologies, Sudbury ore sulfides, and country rocks are plotted as Δ²⁰⁷Pb/²⁰⁴Pb values¹ in Appendix Figure A5 and summarized in the Appendix. Dickin et al. (1992, 1996, 1999) attributed large-scale variations in the Pb isotope compositions of the ores on the North and South Ranges to heterogeneities in the target rocks and small-scale variations to incomplete homogenization of the melt sheet. Darling et al. (2010a, b) attributed vertical variations to variable proportions of melt from different target lithologies, but lateral variations (generally <5 km) to digestion of entrained clasts, fallback material, and local target rocks. O’Sullivan et al. (2016), Kenny et al. (2017), and McNamara et al. (2017) suggested that the melt sheet may have been isotopically more homogeneous and that volatilization of significant (up to 98%) Pb during impact made the system more sensitive to incorporation of Pb during subsequent postimpact thermomechanical erosion of isotopically heterogeneous footwall rocks.

No sublayer inclusions have yet been analyzed for Lu‐Hf isotopes, but available data for offset dikes, main mass lithologies, and exposed and inferred country rocks are plotted as ε_{Hf(1850 Ma)} values in Appendix Figure A6 and summarized in the Appendix. The extremely small variations in the main mass and offset dikes require the melt sheet to have been well homogenized during impact mixing and/or subsequent convection. Kenny et al. (2017) attributed the greater homogeneity of the Hf isotope system to it being much less volatile than Pb and, therefore, much less susceptible to modification during postimpact thermomechanical erosion.

The mass‐dependent ³²S and ³⁴S isotope compositions of main mass lithologies, the sublayer, Sudbury sulfide ores, and overlying and underlying rocks are plotted as δ³⁴S values in Appendix Figure A7 and summarized in the Appendix. There are significant heterogeneities within and between units/zones, but there does not appear to be any correlation between S isotope composition and the overall grade or tonnage of sulfide mineralization: high‐grade and low‐grade deposits have similar ranges in δ³⁴S (Ripley et al., 2015). Sulfur is as volatile as Pb (e.g., Lodders et al., 2009), so it also should have been volatilized, requiring much of the S to be derived from S-bearing footwall rocks (Huronian rocks, East Bull Lake Intrusive Suite intrusions, and Nipissing intrusions) during postimpact thermomechanical erosion. The S contents of marginal quartz diorite vary continuously down to below lower limits of detection (<50 ppm: Wallbridge Mining Company, unpub. data, 2021), consistent with significant S loss.

The mechanics of impact cratering are fairly well understood (see reviews by Grieve, 1987; Melosh, 1989; French, 1998; Osinski and Pierazzo, 2013) and provide important constraints on the formation and localization of the sublayer, footwall breccia, inclusion quartz diorite, and associated Ni-Cu-PGE mineralization. The key constraints are:

1. Impacts involve several discrete stages of differing durations: contact and compression (seconds), volatilization and excavation (seconds to minutes), isostatic rebound and formation of a central uplift (minutes), collapse of the central uplift and crater walls (minutes), crystallization of the impact melt and additional structural readjustments (tens of thousands of years), and tectonic modification and erosion (up to millions of years).

2. Large impacts like Sudbury produce larger proportions of impact melt and volatiles and smaller proportions of impact debris (e.g., Grieve and Cintala, 1997; Cintala and Grieve, 1998). However, abundant debris would have been generated during rebound and collapse.

3. The largest impacts (e.g., Chicxulub, Sudbury, Vredefort) produce impact melt dikes that are tens of meters wide and extend tens of kilometers into underlying rocks. They are interpreted to be emplaced during excavation (e.g., Grant and Bite, 1984; Lightfoot and Farrow, 2002; Murphy and Spray, 2002) and/or rebound (e.g., Wood and Spray, 1998; Tuchscherer and Spray, 2002; Dressler and Reimold, 2004).

4. Thus far, only Sudbury contains magmatic Ni-Cu-PGE mineralization.

Everything that followed impact—including the formation of the sublayer, footwall breccia, and inclusion quartz diorite; ore-localizing embayments, troughs, and funnels; and associated mineralization—depended on the degree of superheat and the amount of impact debris that was produced. Some geologic models suggest that original impact topography was preserved (e.g., Morrison, 1984), implying that the impact melt was never significantly superheated, but most thermal models imply complete digestion of impact debris (e.g., Ivanov and Deutsch, 1999; Prevec and Cawthorn, 2002). Some petrological models suggest slower digestion of impact debris, providing time for gravitational segregation into upper felsic and a lower mafic layers (e.g., Golightly, 1994; Zieg and Marsh, 2005), but the occurrence of melanoritic rocks throughout much of the stratigraphy of the Sudbury Igneous Complex, which was interpreted to likely be fragments of a roof melanoritic sequence that grew from the top and was later disrupted, is inconsistent with such models (Latypov et al., 2019).

Any model for the cooling history of the impact melt must be consistent with the following:

1. Temperatures adjacent to the impactor likely reached 3,700° to 4,700°C (Grieve et al., 1977; Wünnemann et al., 2008), but the majority of the impact melt is interpreted to have been superheated to ~1,700°C (e.g., Ivanov and Deutsch, 1999; Abramov and Kring, 2004; Ivanov, 2005).

2. There are significant amounts of exotic 2.6 to 2.5 Ga zircons in the Onaping Formation (Petrus et al., 2016), but no exotic zircons in quartz diorite, suggesting that they melted at temperatures greater than ~1,670°C (Ostermann et al., 1996).

3. The aphanitic and spherulitic textures of thin (<20 cm) distal (10–15 km) quartz diorite apophyses (e.g., Trill, Hess) are consistent with the impact melt having been superheated when it was emplaced.

4. The sparse inclusions in marginal quartz diorite are locally derived country rocks, not impact debris, and do not include any exotic ultramafic xenoliths, suggesting that all of the impact debris was digested.

5. The abundant inclusions in the protosublayer and proto-inclusion quartz diorite must have been produced after emplacement of marginal quartz diorite, most likely during generation and collapse of the central uplift and crater walls. This would have rapidly cooled the impact melt, but in contrast to the undifferentiated 230-m-thick Manicouagan melt sheet, which has a glassy matrix with ~25% clasts, indicating rapid cooling to the liquidus (Onorato and Uhlmann, 1978), the Sudbury impact melt would have remained above the liquidus (arbitrarily assumed here to be ~1,300°C) during injection of protoinclusion quartz diorite.

6. The usual restriction of inclusion quartz diorite to the cores of offset dikes suggests that it intruded before the quartz diorite margins had completely solidified, which we estimate to have taken ~1 to 5 days (see App. section 6.2). This time is within the residence times of >8-cm-sized ultramafic inclusions and >10-cm-sized felsic or hydrated basaltic inclusions at 1,300°C (see “Convective settling” section below; App. section 5).

7. The exotic ultramafic inclusions in the sublayer and inclusion quartz diorite are rounded. They generally do not exhibit any textural or geochemical evidence of partial melting (Wang et al., 2020), but may have been partly dissolved (see Kerr, 1994a, b).

8. The abundant locally derived inclusions in sublayer and inclusion quartz diorite require significant amounts of thermomechanical erosion and therefore superheat.

9. Systematic gradations in fragment shapes and orientations across the footwall breccia and sublayer indicate that they represent progressive thermomechanical erosion of coherent footwall rock (Gregory, 2005).

10. The exotic ultramafic inclusions in the footwall breccia must have been acquired during interaction with the sublayer.

Together, these constraints suggest that the impact melt was superheated to 1,700°C or more, that it was inclusion-free prior to emplacement of marginal quartz diorite, that it was rapidly cooled by rebound/collapse debris prior to emplacement of inclusion quartz diorite, and that it retained enough heat to thermomechanically erode local footwall rocks to form footwall breccia and add components to the sublayer.

The lower contact of the Sudbury Igneous Complex presently exhibits a significant variety of topography (amplitudes up to 400 m over wavelengths of hundreds of m to a few km and amplitudes up to 1,500 m over a wavelength of 25 km: see reviews by Dreuse et al., 2010; Lightfoot, 2016), most of which has been attributed to primary crater floor topography (e.g., Morrison, 1984; Golightly, 1994; Dreuse et al., 2010; Lightfoot, 2016). However, if the Sudbury Igneous Complex remained superheated, there would be significant amounts of postrebound/collapse thermomechanical erosion, decreasing the height of constructive topographic features (e.g., peak rings: Prevec and Cawthorn, 2002) and increasing the depths of destructive topographic features (e.g., embayments, troughs, and funnels).

Even if superheated by only tens of degrees, the 2- to 5-km-thick impact melt sheet would have convected vigorously (see App. section 5), removing any initial chilled margin at the base of the Sudbury Igneous Complex (see Huppert and Sparks, 1985), melting large amounts of uplift/collapse debris (see below), and melting up to 500 m of the footwall rocks (Prevec and Cawthorn, 2002). The rate of thermal erosion would have varied with time- and location-dependent variations in (1) the thickness, composition, density, viscosity, thermal properties, and convective regime of the impact melt, (2) the geometry of residual and generated footwall topography, and (3) the compositions, densities, thermal properties, and textures of the underlying and overlying rocks. Thermal erosion would have been aided by mechanical erosion of rocks fractured during impact, particularly beneath or adjacent to embayments where dense molten sulfides percolated into fractured footwall rocks, generating the footwall breccia (e.g., Gregory, 2005).

The wide variations in all of these many parameters make it difficult to develop a mathematical model of the erosion history, but we can estimate minimum thermal melting rates for wall rocks and melting times and residence times of variably sized inclusions using the methods of Shibano et al. (2013) and Robertson et al. (2015) (see App. section 5). Thermal erosion rates at 1,700°C (during injection of quartz diorite) and at 1,300°C (estimated during injection of inclusion quartz diorite) are ~0.4 and ~0.04 m d^–1, respectively, for felsic or hydrated mafic rocks, ~0.2 and ~0.006 m d^–1 for dry mafic rocks, and ~0.04 and ~0.0004 m d^–1 for ultramafic rocks (Fig. 3A). At the same two temperatures, small (0.01–0.1 m) felsic and hydrated mafic inclusions will melt in 0.2–1.7 and 4.0–30 min, small mafic inclusions in 0.3–5.4 and 5.8–95 min, and small ultramafic inclusions in 0.3–17 and 5.7–297 min (Fig. 3B). Large (1–10m) felsic and hydrated mafic inclusions will melt in 0.05–0.4 and 0.8–6.6 days, large mafic inclusions in 0.07–1.2 and 1.3–21 days, and ultramafic inclusions in 0.07–3.7 and 1.3–60 days (Fig. 3B). If the inclusions were heavily fractured, they would have been more easily disaggregated, reducing their residence times. Similarly, if the inclusions contained any Fe ± Ni ± Cu sulfides, which would have started melting at 700° to 900°C (depending on composition: e.g., Naldrett, 1969; Tsujimura and Kitakaze, 2004; Helmy et al., 2021), this would also have reduced their residence times.

Mechanical erosion by very dense (~4.2 g cm^–3), very low-viscosity (~0.05 Pa s) molten Fe-Ni-Cu sulfide liquids would have continued until they solidified at ~700° to 800°C (Tsujimura and Kitakaze, 2004; Helmy et al., 2021). Flow rates for molten sulfides into footwall fractures should be much higher than thermal erosion rates of floor rocks by impact melt but will be limited by the ability of the sulfides to enter fractures without freezing. However, the extent of infiltration will increase with time as the rate of thermal erosion decreases and heat is conducted further and further into the footwall rocks, forming mineralized anatectic footwall breccia and underlying mineralized footwall veins and stockworks (Gregory, 2005).

All of the inclusions will float in dense sulfide liquid (Fig. 3A), indicating that the inclusions that remain were frozen in place as the sulfide solidified.

The settling of sulfide droplets during convective flow and gravity flow of the protosublayer (and protoinclusion quartz diorite) and protoinclusion semimassive sulfides along the lower contact of the melt sheet are evaluated in the Appendix and summarized below.

Estimated settling velocities and residence times of various sized inclusions and sulfide droplets are compared to the mean mid-depth convective velocity of the impact melt in Figure 4. The convective velocity is greater than settling velocities of <5-cm sulfide droplets and <80-cm inclusions at 1,300°C. However, inclusions and sulfide droplets will still settle out of even a vigorously convecting magma because convective velocities approach zero at the bases of convection cells (see Martin and Nokes, 1989). Small sulfide droplets (20–50 μm: the sizes expected if they exsolve from the impact melt during cooling) will also settle at the bottom of convection cells but at much slower rates (Fig. 4). Importantly, the estimated 1- to 5-day interval between emplacement of marginal quartz diorite and inclusion quartz diorite is much shorter than the 30–200 years required for 20- to 50-μm sulfide droplets to exsolve and settle through 1 km of norite at 1,300°C (see Fig. 4; App. Table A3), as required in the sulfide exsolution model.

The estimated residence times of small inclusions (1–10 cm), large inclusions (1–10 m), and sulfide droplets (0.5–2 cm) are in the range of 0.9 to 312 days, 13 to 45 min, and 122 to 183 days at 1,300°C (Fig. 4B). Although sulfides will remain immiscible, we have shown above that inclusions will be assimilated by the superheated impact melt as they settle. It takes a few minutes to melt small inclusions (before they settled to the bottom), but a few days to melt large inclusions (Fig. 3). The combined effects of differential settling and preferential assimilation of felsic > mafic > ultramafic inclusions at ~1,300°C explains the close association of mafic-ultramafic inclusions and sulfide droplets in sublayer and inclusion quartz diorite.

Gravity flow of mixtures of impact melt, inclusions, and sulfide liquid has been considered in the Appendix (section 5.2). Depth-averaged horizontal flow velocity is strongly dependent on the pore pressure ratio (ratio of pore pressure to ambient pressure = overpressure), which is unknown in our system and inhibits us from making velocity estimates, but it is clear that with all else equal, velocity is directly proportional to increasing reduced bulk density (Δρ = bulk density – impact melt density) and decreasing bulk viscosity. The protosublayer would have been denser but more viscous than overlying impact melt, whereas protoinclusion semimassive sulfides would have been much denser but much less viscous than overlying impact melt (App. Figs. A8, A9). The closest analogues, submarine debris flows, have ratios of reduced density to viscosity ratios (Δρ/η) in the range of 850 to 1,722, which are significantly greater than the protosublayer (Δρ/η ~ 0.02–0.07) but overlap the lower limit of protoinclusion-semimassive sulfide (Δρ/η ~ 1,229–24,300) (Fig. 5; App. Table A4). Because the velocity of gravity flow is proportional to Δρ/η, we can predict that protoinclusion-semimassive sulfide will flow down even shallow topographic gradients into embayments, troughs, and funnels and/or accumulate in hydraulic jumps where there are transitions in slope or topography (Fig. 6), just as mixtures of water and sediment flow down slopes in the sea floor as submarine debris flows. However, the protosublayer is much less likely to form gravity flows except on steep slopes.

The East Bull Lake and Nipissing intrusive bodies contain abundant disseminated (James et al., 2002; Holwell and Keays, 2014) and locally net textured (e.g., Shakespeare: Sproule et al., 2007) Fe‐Cu‐Ni sulfides and must be the ultimate source of the sulfides (e.g., Lightfoot et al., 1997; Keays and Lightfoot, 2004). However, few inclusions (most anteliths, but only rare local xenoliths and no exotic xenoliths) contain any sulfides (see App. section 2), there is no evidence of Ni depletion in any of the olivines in the mafic‐ultramafic inclusions (App. Fig. A2), there is no positive correlation between Ni and S in the majority of mafic‐ultramafic inclusions, except for some sulfide‐phlogopite–bearing olivine melanorite anteliths (App. Fig. A2), and local thermomechanical erosion cannot explain the exotic inclusions. This means that there is no genetic relationship between the exotic ultramafic inclusions and Ni-Cu-PGE mineralization.

Another relationship might be hydraulic equivalence (similar fluid dynamic behavior under the same conditions: e.g., Rittenhouse, 1943). However, the size range of molten sulfide blebs in the sublayer is fairly narrow (<1–2 cm), whereas the size range of inclusions is much greater (up to 10s of m). Both would settle more rapidly than mafic-intermediate-felsic inclusions (Fig. 3A), but only a small size range of ultramafic inclusions would be hydrodynamically equivalent.

The most likely influence is survivability. After sulfide saturation was reached, any excess sulfides would have been immiscible in the impact melt. Olivine-bearing mafic and ultramafic inclusions might dissolve at more or less the same rate as mafic-intermediate felsic inclusions but would melt much more slowly (Fig. 3).

The inferred timing relationships of impact cratering, offset dikes, the sublayer, inclusion quartz diorite, and footwall breccia are summarized in Table 2 and Figure 7. The following temporal sequence of events is proposed for the formation of the Sudbury structure, with emphasis on the timing of generation of the inclusion‐rich sublayer, inclusion quartz diorite, footwall breccia, and associated Ni‐Cu‐PGE sulfide mineralization, shown graphically in Figure 8.

An ~15-km‐diameter comet (e.g., Pope et al., 2004; Petrus et al., 2015) impacted the Sudbury area at 1850 Ma, volatilizing up to 98% of Pb, 78% of Zn, 90% of Cs, 30% of Rb, significant amounts of Cd (O’Sullivan et al., 2016; Kenny et al., 2017; McNamara et al., 2017; Kamber and Schoenberg, 2020), almost certainly also S-Hg‐Tl‐Se‐Sn‐Te‐Bi, and possibly also some Sb-Ag‐Cu‐Au‐As (Lesher, 2019a) relative to less volatile chalcophile (Co‐Ni‐PGE) and lithophile (e.g., Th‐Nb‐Ta‐Hf‐ Zr‐Y‐rare earth) elements, producing a superheated devolatilized impact melt with an initial temperature of ~1,700°C (Ivanov and Deutsch, 1999).

The excavation stage would have produced an upward‐directed pressure gradient that expelled impact melt (preserved 500–900 km to the west: see review by Lightfoot, 2016) and airborne debris (proto-Onaping Formation), leaving inclusion‐free melt in the transient cavity. This is the only time when inclusion‐ and sulfide‐free impact melt was present and could be injected into radial and concentric fractures to form the inclusion‐ and sulfide‐free marginal phases of offset dikes.

The compositions of marginal quartz diorite were modified by contamination during emplacement but are close to the composition of the original impact melt (Lightfoot et al., 1997).

Excavation would have been rapidly followed by uplift of the deeper central parts (source of high‐pressure exotic ultramafic inclusions), formation of a central uplift, and rapid collapse to form one or more rings separated by relatively flat circular basins (e.g., Grieve et al., 1977, 1981; Golightly, 1994; Lightfoot, 2016), albeit likely much more subdued because of the very large size of Sudbury (R. Grieve, pers. commun., 2021). This would have been rapidly followed by equally rapid collapse of the peak rings and crater walls (source of local ultramafic inclusions). Felsic and intermediate inclusions melted, leaving mainly mafic‐ultramafic xenoliths. Importantly, most of the footwall rocks contained minor sulfides (Huronian basalts and sediments), many contained significant amounts of sulfides (e.g., Nipissing Intrusive Suite: Lightfoot and Naldrett, 1996; East Bull Lake Intrusive Suite: James et al., 2002), and some contained Ni-Cu-PGE ores (e.g., Shakespeare: Sproule et al., 2007), leading to rapid sulfide saturation and formation of sulfide xenomelts, which were injected—along with refractory exotic and local mafic‐ultramafic xenoliths—into the still‐hot cores of the quartz diorite dikes, forming the observed nested quartz diorite/inclusion quartz diorite dikes and offset ores. Inclusion quartz diorite was isolated from the impact melt, so it cooled and crystallized rapidly, was able to assimilate few inclusions, and was able to accumulate little orthopyroxene-plagioclase and expel little residual silicate melt.

The impact melt continued to thermomechanically erode local footwall rocks until it reached the liquidus. Where assimilated rocks were more magnesian than the impact melt, the melt reached the liquidus prior to the remainder of the impact melt and crystallized a layer of olivine melanorite that was subsequently disrupted to form anteliths in the Foy and Whistle embayments (e.g., Corfu and Lightfoot, 1996; Prevec et al., 2000; Wang et al., 2020).

Any mixtures of exotic and local inclusions and sulfide xenomelts that were not mobilized into the offset dikes remained along the basal contact and flowed into topographic lows under the influence of convection currents and/or gravity. The sublayer remained in contact with the overlying impact melt and so crystallized relatively slowly, remained in communication with overlying impact melt, was able to assimilate most felsic‐intermediate inclusions and parts of mafic‐ultramafic inclusions, and accumulated orthopyroxene-plagioclase to form a noritic matrix. The typical contact ore segregation profile grading downward from disseminated sulfide in sublayer norite through inclusion-bearing net-textured/semimassive sulfide ores to semimassive sulfide ores (Souch et al., 1969) can be attributed to gravitational segregation (e.g., Naldrett, 1973; Ussleman et al., 1979) and/or downward percolation (e.g., Chung and Mungall, 2009; Barnes et al., 2016).

The footwall breccia appears to have formed over the widest range of temperature and time, but the early products of thermomechanical erosion were probably incorporated into the sublayer, and what is present now represents the last stages of partial melting of footwall rocks. The footwall breccia contains no quartz dioritic or noritic matrix, so it must have acquired local xenoliths by incorporating them from the sublayer and/or by thermomechanical erosion of footwall lithologies, and exotic xenoliths by incorporating them from the sublayer. When the temperature of the impact melt was above the liquidus, the rate of thermomechanical erosion would have exceeded the rate of heat conduction into the footwall rocks (see discussion by Huppert and Sparks, 1985), continuing to generate the sublayer, but as the melt reached the liquidus, thermomechanical erosion ceased, allowing heat to be conducted into the footwall rocks, forming anatectic footwall breccia that locally crosscuts the sublayer and lower main mass norites.

Where assimilated footwall rocks were orthopyroxene rich, the melt crystallized and accumulated orthopyroxene‐chromite then orthopyroxene‐plagioclase, forming North Range mafic norite (longer period of orthopyroxene‐(chromite) accumulation) and South Range quartz‐rich norite (shorter period of orthopyroxene‐(chromite) accumulation). Where the assimilated footwall rocks were mafic to intermediate, the melt crystallized and accumulated cotectic orthopyroxene‐plagioclase, forming North Range felsic norite and South Range norite.

As the main mass continued to crystallize, heat was slowly conducted increasing distances into the footwall rocks, allowing fractionated residual sulfide melts to percolate increasing distances into fractured footwall rocks (Gregory, 2005), leaving behind MSS‐rich cumulates in the sublayer, footwall breccia, and contact ores (e.g., Li et al., 1992; Mungall, 2007).

This paper is dedicated to the memory of Prof. A.J. Naldrett, who was a supervisor and mentor to PCL, a career‐long mentor to CML, and an inspiration to YW and EFP. In 1984 he published a schematic cross section illustrating an intrusive model for the emplacement of the Sudbury Igneous Complex. A very few years later, he embraced the origin of the complex as an impact melt sheet, but he continued to emphasize the importance of understanding the suite of mafic‐ultramafic inclusions associated with the mineralized embayment and offset structures. By the time of his 2004 textbook, he wrote, “In some areas the Sublayer is characterized by inclusions. These can be divided into two groups, those of obviously local derivation, and those composed of mafic and ultramafic rocks, many of which do not outcrop in the immediate Sudbury area” (Naldrett, 2004, p. 422). He accepted that the complex was a product of impact melting, and he advocated for detailed studies of the geochemistry of the Sudbury Igneous Complex. He thought that the origin of these inclusions provided a window into the formation of the melt sheet and the composition of the target rocks. His ideas were pivotal at the inception of the project giving rise to this paper. Financial support was provided by a Natural Sciences and Engineering Research Council of Canada Discovery grant (203171‐2012) to CML, a China Scholarship Council award to YW, and a Canada First Research Excellence Fund grant to CML and others. Vale Canada Ltd. kindly provided access to drill core samples and their thin section library. We are grateful to Lisa Gibson, Noëlle Shriver, and other Vale geologic staff for assistance with some of the sampling, to Richard Grieve for discussions on impact cratering, and to Rais Latypov, Steve Barnes, and an anonymous Economic Geology reviewer for insightful and helpful reviews on the manuscript, which led to many improvements. We thank Dr. Eduardo Mansur for handling of our manuscript and for patiently waiting for our revision. This is Mineral Exploration Research Centre Metal Earth contribution MERC‐ME‐2021‐42.

Yujian Wang received her Ph.D. degree at Laurentian University in 2019. Her Ph.D. project was related to the petrology, geochemistry, and petrogenesis of mafic-ultramafic inclusions in the sublayer and offset dikes and the relationship with the massive Ni-Cu-PGE sulfide mineralization in the Sudbury Igneous Complex, Canada. She has published the relevant work in Geology, Journal of Petrology, and EconomicGeology. Yujian is now an associate professor at China University of Geosciences (Beijing). She is currently working on the mantle petrology and geochemistry, radiogenic isotope geochemistry, and geodynamics of convergent plate tectonics.