Synchronous Emplacement of the Anorthosite Xenolith‐Bearing Beaver River Diabase and One of the Largest Lava Flows on Earth

New geochronologic and paleomagnetic data from the North American Midcontinent Rift (MCR) reveal the synchronous emplacement of the Beaver River diabase, the anorthosite xenoliths within it, and the Greenstone Flow—one of the largest lava flows on Earth. A U‐Pb zircon date of 1091.83 ± 0.21 Ma (2 σ ) from one of the anorthosite xenoliths is consistent with the anorthosite cumulate forming as part of the MCR and provides a maximum age constraint for the Beaver River diabase. Paired with the minimum age constraint of a cross‐cutting Silver Bay intrusion (1091.61 ± 0.14 Ma; 2 σ ), these data tightly bracket the age of the Beaver River diabase to be 1091.7 ± 0.2 Ma (95% CI), coeval with the eruption of the Greenstone Flow (1091.59 ± 0.27 Ma; 2 σ )—which is further supported by indistinguishable tilt‐corrected paleomagnetic pole positions. Geochronological, paleomagnetic, mineralogical and geochemical data are consistent with a hypothesis that the Beaver River diabase was the feeder system for the Greenstone Flow. The large areal extent of the intrusives and large estimated volume of the volcanics suggest that they represent a rapid and voluminous ca. 1,092 Ma magmatic pulse near the end of the main stage of MCR magmatism.

. (a) Geologic map of exposures of Midcontinent Rift volcanics and intrusives in the western Lake Superior region. The Greenstone Flow (purple) of the Portage Lake Volcanics (red) outcrops throughout the Keweenaw Peninsula and Isle Royale. (b) Regional map of paleomagnetic and geochronologic sites in the southern Beaver Bay Complex (south BBC). Note that paleomagnetic site AX16 and geochronology sample MS99033 are from the same anorthosite xenolith. The geochronology sample numbers in (a) and (b) correspond to those in Figure 4. (c) Regional map of paleomagnetic sites in the eastern Beaver Bay Complex (east BBC). The xenolith at Carlton Peak is  E 100 m in diameter. The younger Schroeder-Lutsen basalt of the North Shore Volcanic Group (NSVG) is lying unconformably atop the Beaver River diabase and other NSVG units. The nomenclature of the "southern" and "eastern" Beaver Bay Complex follows Miller and Chandler (1997). FTMD-Finland tectonomagmatic discontinuity, traced out by the dashed black line. Bedrock geology is from Miller et al. (2001) and Jirsa et al. (2011).
(light green bar), the Greenstone Flow from the 1.1 Ga Midcontinent Rift stands amongst the giant lava flows from the Deccan Traps and Columbia River basalts. M-R Traps = Mahabaleshwar-Rajahmundry lava flow in the Deccan Traps. Volume estimates from Self et al. (2008), Bryan et al. (2010), Longo (1984), and Doyle (2016). 10.1029/2021GC009909 5 of 22 recorded by MCR rocks that were emplaced a few million years apart (Swanson-Hysell et al., 2019). In this study, we present paleomagnetic data from the anorthosite xenoliths and the host Beaver River diabase. Data from the xenoliths give equivalent directions to the host diabase (Figures 7 and 8), indicating that they were heated above the Curie temperature of magnetite and acquired a thermal remanent magnetization when they cooled within the diabase. This thermal history is consistent with thermal diffusion modeling of the xenoliths (Figure 9). The paleomagnetic data can be compared to data from the Greenstone Flow to further test the hypothesis that they are synchronous. The resulting paleomagnetic pole positions can also be compared to the synthesized Laurentia APWP to obtain chronological constraints ( Figure 8).
Here, by integrating the geochronologic and paleomagnetic perspectives with previous lithologic and geochemical analyses (Doyle, 2016;Miller & Chandler, 1997), we show that these data are consistent with the Beaver River diabase network acting as the feeder system for the Greenstone Flow of the PLV. Alternatively, they could both be the distinct manifestations of magmatism from a similar source. Regardless, their shared geochemical signatures and the inference of giant magma conduits that transported large anorthosite xenoliths characterize a period of ca. 1,092 Ma voluminous magmatic activity (based on 206 Pb/ 238 U zircon dates; Figure 1).

Beaver Bay Complex and Related Rocks of NE Minnesota
The North American MCR is a failed intracontinental rift where protracted magmatic activity lasted from ca. 1,109 to ca. 1,084 Ma (Swanson-Hysell et al., 2019). MCR rocks extensively outcrop in today's Lake Superior region, with the total extent traceable by arcuate magnetic and gravity anomalies that extend to the southwest to Kansas, and to the southeast, to southern Michigan (Hinze & Chandler, 2020). Previous studies have divided magmatic activity in the rift into four stages based on interpreted changes in relative magmatic volume and the nature of magmatism: early (  E 1,109-1,104 Ma), latent (  E 1,104-1,098 Ma), main (  E 1,098-1,090 Ma), and late (  E 1,090-1,083 Ma) (Miller & Nicholson, 2013;Vervoort et al., 2007). In northeastern Minnesota, the Early Gabbro Series and the Felsic Series rocks of the Duluth Complex and reversed-polarity lavas of the lower North Shore Volcanic Group (NSVG) were emplaced during the early stage. The more voluminous Duluth Complex Layered Series and the plagioclase-rich Anorthositic Series, together with an associated  E 8 km thick extrusive volcanic sequence of the NSVG, were rapidly emplaced about 10 myr later at ca. 1,096 Ma during the main stage (Paces & Miller, 1993;Swanson-Hysell et al., 2021).
The Beaver Bay Complex, which sits stratigraphically above the Duluth Complex, is another intrusive complex that resulted from main stage magmatism. The exposed area of the Beaver Bay Complex is  E 1,000 km 2 where it has been mapped along the northwestern shore of Lake Superior in northeastern Minnesota (Figure 1). The Beaver Bay Complex is a multi-phase, composite intrusive complex that intrudes parts of the NSVG (Figure 1; Miller & Chandler, 1997;Swanson-Hysell et al., 2021). Distinct from the deep plutonic intrusions of the Duluth Complex, the majority of the Beaver Bay Complex is formed of hypabyssal intrusions that were emplaced as dikes and sills at shallow depths (Miller & Chandler, 1997). Most of the Beaver Bay Complex intrusions are dioritic to gabbroic in composition (Miller & Chandler, 1997). The main lithology of the Beaver River diabase dikes and sills network within the Beaver Bay Complex is an ophitic olivine gabbro ( Figure 2), but in wider areas of dikes and the upper parts of thick sills, this rock type can abruptly transition into intergranular olivine oxide gabbro, then into subprismatic (and commonly foliated) ferrogabbro, and finally into granophyric monzodiorite. The more evolved and later emplaced components of the Beaver River diabase network are commonly distinguished as the Silver Bay intrusions in the southern Beaver Bay Complex (Figure 1). Overall being intermediate in composition, the Silver Bay intrusions lithologies range from ophitic olivine gabbro to ferrogranite (Shank, 1989). Field mapping by Miller et al. (1994) found intrusive relationships between the Silver Bay intrusions and the Beaver River diabase. Angular inclusions of the host Beaver River diabase within marginal zones of the Silver Bay intrusions led Miller and Chandler (1997) to interpret that the Silver Bay intrusions intruded after the diabase crystallized.
One distinctive feature of the Beaver River diabase is its inclusions of anorthosite xenoliths. In the southern part of the Beaver Bay Complex, the Beaver River diabase occurs as dikes and sills, typically including anorthosites with various sizes ranging from centimeters to over 150 m (Figures 1 and 2; Grout & Schwartz, 1939;Morrison et al., 1983). The diabase in this region intrudes the Palisade rhyolite of the NSVG (Figure 1), which has a 206 Pb/ 238 U date of 1093.94  E 0.28 Ma (2  E analytical uncertainty is presented for CA-ID-TIMS dates throughout this work; Swanson-Hysell et al., 2019). The Beaver River diabase is locally intruded by the Silver Bay intrusions (Figure 1). An aplite unit within the granophyre zone of one of these Silver Bay intrusions has a 206 Pb/ 238 U date of 1091.61  E 0.14 Ma (Swanson-Hysell et al., 2019). Another arcuate, sill-like diabase body mapped as the Beaver River diabase outcrops along the eastern part of the complex (Figure 1; Miller & Chandler, 1997). The diabase composition there is similar to that in the south and it also contains large anorthosite xenoliths with dimensions that exceed 100 m at Carlton Peak ( Figure 1). The Beaver River diabase in the northern part of the complex, near the Houghtaling Creek area, typically forms narrow, near-vertical dikes instead of sheets in the southern and eastern regions (Figure 1; Miller et al., 1994). The diabase in this region only locally contains xenoliths of anorthosite.
Hundreds of anorthosite xenoliths have been recognized and mapped within the Beaver River diabase (Figure 1). Many hill tops in the Beaver Bay Complex, such as Carlton Peak and Britton Peak, are large anorthosite blocks (which lead Lawson, 1893 to erroneously conclude that they were relict Archean topography). Later work established the anorthosite blocks as xenoliths, which are now extensively documented through geologic mapping of the region (Figure 1; Boerboom, 2004;Boerboom et al., , 2007Miller, 1988;Miller & Boerboom, 1989;Miller et al., 2001) and outcrop-scale exposures ( Figure 2). In the field, the anorthosites typically appear as subrounded to rounded, light-colored, translucent blocks that are in sharp contact with the hosting diabase ( Figure 2). They also occur as exposures whose contact with the diabase is covered (Figure 2). Grout and Schwartz (1939) suggested that the rounded anorthosites are the result of abrasion during transportation as they were entrained by the diabase. While the Beaver River diabase is chilled against the NSVG lithologies that it intrudes, the diabase is not chilled against the margin of the anorthosite xenoliths (Miller & Chandler, 1997;Morrison et al., 1983). The lack of chilled contacts is consistent with the anorthosite being at elevated temperatures and cooling at the same time as the diabase magma ( Figure 9).
The anorthosite xenoliths are dominantly monomineralic plagioclase that has an average anorthite content of  E 70% (Doyle, 2016;Morrison et al., 1983). Interstitial pyroxene and olivine are present in minor concentrations in the xenoliths. Within the Carlton Peak anorthosite xenolith, up to 10 cm oikocrysts of olivine and pyroxene can occur. Nevertheless, the overall olivine content in the anorthosites is low. Interstitial titanomagnetite-ilmenite intergrowths that exceed 100  E m can be found with microscopy and  E 20  E m Fe-Ti oxide grains can be detected with scanning electron microscopy ( Figure 2). Based on textural differences, Morrison et al. (1983) divided the anorthosite xenoliths into four groups: one group which typically have well-developed granoblastic texture characterized by equigranular plagioclase crystals; another group which have interlocking, lath-shaped plagioclase crystals; an intermediate group which can have both granoblastic texture and interlocking plagioclase laths; and a brecciated group that have brittle deformation textures superposed on pre-existing textures.

Portage Lake Volcanics and the Greenstone Flow
The Portage Lake Volcanics is a  E 5 km thick, normally magnetized, dominantly olivine basalt to andesite volcanic succession that outcrops in northern Michigan (particularly along the Keweenaw Peninsula) as well as on Isle Royale (Figure 1; Cannon & Nicholson, 2001;Green, 1982;Huber, 1973). The Greenstone Flow of the Portage Lake Volcanic Group has been recognized as one of the largest lava flows on earth (Figures 1 and 3). It outcrops as the main ridge along the Keweenaw Peninsula and Isle Royale (Figure 1). The flow can be correlated between the two outcrop regions on the basis of geochemical, petrographic, and paleomagnetic similarity of the flow itself and the flows above and below (Longo, 1984). In both outcrop regions, the Greenstone Flow is underlain by conglomerate and overlain by pyroclastic breccia (Huber, 1973;Lane, 1911). On the Keweenaw Peninsula, the Greenstone Flow is exposed over 90 km with a range of thickness from  E 100 m to a maximum thickness of over 450 m, dipping to the northwest (Figure 1; White, 1960). On Isle Royale, the Greenstone Flow has a range of thickness from  E 30 m to a maximum thickness of about 250 m, dipping toward the southeast (Figure 1; Huber, 1973). More recently, Doyle (2016) estimated that the total aerial extent of the Greenstone Flow could be up to  E 20,000 2 km E by connecting it to the region of the Beaver Bay Complex. Taking thickness range of 100-300 m, Doyle (2016) estimated a total volume of According to the mineralogical and textural attributes, the Greenstone Flow can be divided into four zones from bottom to top-a lower ophitic zone, a "pegmatoid" or heterolithic zone, an upper ophitic zone, and an amygdaloidal zone (Cornwall, 1951). The heterolithic zone contains lenses to layers of coarse-grained granophyric gabbro that are referred to in the literature as "pegmatoid." Zircon crystallized in these layers have enabled the heterolithic zone to be targeted for U-Pb geochronology (Davis & Paces, 1990;Swanson-Hysell et al., 2019). A 206 Pb/ 238 U zircon date of 1091.59  E 0.27 Ma for the Greenstone Flow was developed from a sample from this zone in Swanson-Hysell et al. (2019). The Greenstone Flow is typically interpreted to represent emplacement of a single body of magma that then underwent in situ differentiation (Davis & Paces, 1990;Huber, 1973). Doyle (2016) favored a distinct model in which the Greenstone Flow is a composite unit, which they interpret to be indicated by lithologic zonation of ophitic basalt forming the upper and lower zones and an interior zone composed of prismatic ferrogabbro to granophyric monzodiorite. They envision emplacement of the Greenstone Flow started with voluminous eruption of olivine tholeiitic magma, forming the ophitic zones which while still crystallizing further inflated due to subsequent injection of a more evolved basaltic magma to form intergranular gabbro in the heterolithic zone. They considered this progression to be more consistent with observed abrupt lithologic changes from the ophitic zone to the heterolithic zone over centimeter to meter scales, inclusion relationships between evolved and ophitic Greenstone Flow lithologies, and remnant blocks of initially crystallized ophitic basalt interlayered with evolved lithologies within the heterolithic zone which contains the pegmatoids. In both the Doyle (2016) model of multiple magma injections and the earlier models of in situ differentiation, it is migration of the most evolved and volatile-rich melts within the interior of the flow in the final stages of flow crystallization that led to the formation of some aplite dikes and the coarsest segregations containing granophyre. Both models also invoke a single basaltic parental magma with distinction of where differentiation occurred either through fractional crystallization in an evolving magma chamber or solely within a single, very thick flow.

Zircon Geochronology and Geochemistry
A sample of an anorthosite xenolith within the Beaver River diabase was collected for U-Pb geochronology along Hwy 61 across from the Silver Bay taconite plant (MS99033; 91.26358°W 47.28888°N; Figure 1). This sample comes from the same xenolith sampled for paleomagnetic study as site AX16 which has an exposed diameter of 27.5 m ( Figure 2). Thin sections were made from the geochronology sample as well as multiple paleomagnetic cores. As is shown in Figure 2f, plagioclase in this anorthosite xenolith have both equigranular crystals displaying a granoblastic texture and lath-shaped crystals displaying an interlocking texture. The occurrence of both textures is consistent with an interpretation that this anorthosite xenolith formed under elevated temperatures and experienced heating after initial crystallization.
Zircons were separated from a kilogram of the anorthosite using common mineral separation methods (Supporting Information S1). The separated zircons were subhedral to anhedral crystals (z1-z4) and platy fragments (z5-z8). The subhedral to anhedral crystals are consistent with intercumulus crystallization within an adcumulate with platy fragments being a common zircon morphology within anorthosites (e.g., sample AS3 of the Duluth Complex anorthositic series; Schmitz et al., 2003). Eight chemically abraded zircons were analyzed by isotope dilution-thermal ionization mass spectrometry (ID-TIMS) in the Boise State Isotope Geology Laboratory using EARTHTIME tracer solutions (Condon et al., 2015; Supporting Information S1). Both zircon morphologies yield indistinguishable dates. Using six of these single-grain dates (and excluding two due to interpreted Pb-loss) results in a weighted mean 206 Pb/ 238 U date of 1091.83  E 0.21/0.37/1.15 Ma (analytical/analytical + tracer/analytical + tracer + decay uncertainty; Figure 4).
This date provides a tight constraint on the age of the Beaver River diabase. Previously, the maximum age constraint for the Beaver River diabase came from the relationship that it cross-cuts the Palisade rhyolite of the NSVG which has a 206 (Fairchild et al., 2017), cross-cut the Beaver River diabase. These dates constrain the diabase to have been emplaced between 1091.83  E 0.21 and 1091.61  E 0.14 Ma (Figure 4). Assuming a uniform probability of diabase emplacement between the anorthosite and aplite dates and their normal distributed uncertainties, a 95% confidence interval on the age of the diabase can be estimated by Monte Carlo simulation. This analysis gives an age for the diabase of 1091.7  E 0.2 Ma (95% CI).
An additional 15 zircons were characterized using cathodoluminescence (CL) imaging and laser ablation-inductively coupled plasma mass spectrometry (LA-ICPMS), with methods and instrumentation described in Supporting Information S1. CL images reveal internal planar zones of variable brightness, often with darker interior zones and outer brighter zones ( Figure 5). All crystals exhibit sharp, micron-scale transitions between zones, and LA-ICPMS analyses quantify CL brightness as correlated with rare earth elements (REEs) content. REE patterns in zircons exhibit a significant chondrite-normalized negative Eu anomaly ( Figure 6). The Ti-in-zircon thermometer gives a range of estimated zircon crystallization temperatures from 998°C to 860°C with a mean of  E 950°C (Ferry & Watson, 2007; Supporting Information S1). Decreasing temperatures are correlated with deepening of the negative Eu anomaly and increasing incompatible trace element (e.g., Hf, Th) incorporation into zircon. These data are consistent with a model of magmatic zircon crystallizing from cooling and fractionating interstitial residual melt within the cumulate plagioclase framework.

Paleomagnetism
We collected paleomagnetic cores that are 2.5 cm in diameter along the southern and eastern Beaver Bay Complex with a particular focus on acquiring paired sites of anorthosite xenoliths and their local diabase hosts. Sample cores were collected using a hand-held gasoline-powered drill and were oriented using a magnetic compass as well as a sun compass when possible. Sun compass orientations were preferentially used for determining the sample azimuth. Typically, 7-10 cores were collected for each anorthosite xenolith and their diabase hosts. A total of 17 diabase and 22 anorthosite sites were sampled (Table 1). A table that summarizes the measured dimensions of each anorthosite xenolith sampled and the distance between each anorthosite paleomagnetic site and closest diabase host site is provided in Supporting Information S1.
Samples underwent step-wise demagnetization and analyses in the magnetically shielded room at the UC Berkeley Paleomagnetism Lab. Seven sites from the Beaver River diabase underwent alternating field (AF)  (Swanson-Hysell et al., 2019, 2021. These high-precision dates are consistent with field observations that the Beaver River diabase cross-cuts the Palisade rhyolite (dark green) and is cut by the Silver Bay intrusions (pink). The estimated age of the Beaver River diabase from these constraints is shown by an orange box representing the 95 % E confidence interval. Each vertical bar corresponds to one 206 Pb/ 238 U date from a single zircon crystal. The translucent bars represents zircons with interpreted Pb loss that are not included in the weighted mean age calculations. Horizontal lines and gray boxes represent weighted mean 206 Pb/ 238 U dates and their analytical uncertainty. The numbers of each geochronology sample correspond to those in Figure 1 where locations of these samples are shown.
4 5 3 2 1 demagnetization with peak fields from 1 to 130 mT. An ASC TD-48SC thermal demagnetizer was used to demagnetize 10 diabase sites and all 22 anorthosite sites in a step-wise manner, with reduced step increments between 540°C and 585°C. The typical magnetic field inside the shielded room is  E 500 nT and the field inside the thermal demagnetizer chamber is  E 10 nT. The quartz glass sample rod of the UC Berkeley system is typically measured at 5   12 2 10 Am E . All remanence measurements were made on a 2G Enterprises DC-SQUID superconducting rock magnetometer equipped with inline AF coils and an automated sample m. Note that grain 1 (corresponding to spot 1) has a platy morphology, while the rest of the grains are subhedral to anhedral.

Figure 6.
Rare earth element (REE) analyses for four plagioclase crystals and for 15 zircons from geochronology sample MS99033 (anorthosite xenolith site AX16) developed by inductively coupled plasma mass spectrometry. All data are chrondrite-normalized (Sun & McDonough, 1989). All data can be accessed in Supporting Information S1 and Table S1.  Figure 7. Example orthogonal vector demagnetization diagrams for diabase and anorthosite specimens. Anorthosite site AX1 is a xenolith within the diabase sampled as BD1. Similarly, AX22 is from a xenolith within the BD10 diabase. Both AF and thermal demagnetization show dominantly univectoral decay of characteristic remanent magnetizations (ChRM) toward the origin after removal of minimal secondary components. The data show very similar ChRM directions between the paired diabase and anorthosite xenoliths sites. Representative magnetization intensity versus thermal demagnetization step plots are paired with orthogonal vector plots for specimen BD10-3a and AX22-3a. changer system. The PmagPy software package was used to implement least-square fits to specimen demagnetization data (Tauxe et al., 2016). Measurement level data are available within the MagIC database (https://earthref.org/MagIC/doi/10.1029/2021GC009909) For both the diabase and anorthosite demagnetization, principal component analyses show that an origin-trending characteristic remanent magnetization (ChRM) can be isolated after the removal of a minimal secondary component during the first few low coercivity (  E 10 mT) or low temperature (  E 200°C) demagnetization steps (Figure 7). The ChRMs typically unblock through thermal demagnetization steps from  E 500°C to  E 580°C, consistent with the component being held by low-titanium titanomagnetite. We interpret this component as a primary remanent magnetization acquired during the emplacement and cooling of the Beaver River diabase.
The site mean paleomagnetic directions are shown in Table 1. We present both AF and thermal demagnetization results for the Beaver River diabase as both methods are effective in removing the secondary components and isolating the coherent and univectoral ChRM. Based on specimen-and site-level demagnetization behavior and the proximity between paired paleomagnetic sites of the anorthosite xenoliths and the diabase, we grouped the anorthosite xenoliths and their diabase hosts into individual cooling units and calculated a paleomagnetic pole position from the mean of the cooling unit virtual geomagnetic poles (Figure 8). Tilt-correcting the paleomagnetic directions to paleohorizontal is necessary for developing accurate paleomagnetic poles from the diabase and the anorthosite xenoliths to be compared to the Keweenawan Track APWP (Figure 8, Swanson-Hysell et al., 2019). For intrusive igneous rocks, tilt corrections can be difficult to constrain due to the lack of a clear paleohorizontal reference. Many paleomagnetic studies of MCR intrusive rocks in the Lake Superior region did not apply tilt corrections to their data (e.g., Beck, 1970;Beck & Lindsley, 1969;Books et al., 1966). However, we can determine the structural orientation of the Beaver River diabase using the abundant igneous fabric orientations measured on the diabase as well as bedding orientations measured from adjacent volcanic units (Boerboom, 2004;Boerboom et al., , 2007Miller et al., 2001). We compile the igneous layering measurements from the Beaver River diabase and the volcanic bedding orientations from the Schroeder-Lutsen basalt which is overlying the Beaver Bay Complex (Figure 1). Despite the uncertainties associated with using igneous fabrics orientations as paleohorizontal references, the mean tilt orientations of the fabrics of the Beaver River diabase and the volcanic bedding orientations of the Schroeder-Lutsen basalt are similar (diabase overall dip directiondip: 128.5-10.2; basalt dip direction-dip: 142.2-13.6). We combine the structural measurements from the Beaver River diabase and the Schroeder-Lutsen basalt and derived two sets of tilt corrections for the paleomagnetic directions of the diabase and anorthosite (dip direction-dip in the southern Beaver Bay complex: 128.7-12.9; in the eastern Beaver Bay Complex: 145.6-13.1, Supporting Information S1). The advantage of using the structural orientations from the Schroeder-Lutsen basalt is that the arcuate shape of the Beaver River diabase intrusions is nicely captured by the variation of lava dip directions while the dip angles of the basalt and diabase are very similar (Figure 1).
The tilt-corrected ChRMs in both lithologies are west-northwest and down, yielding good specimen-and site-level consistency (Figures 7 and 8). Close directional similarities between each anorthosite xenolith and their host diabase are supported by 9 out of a total of 17 diabase-anorthosite paleomagnetic site pairs passing a common mean test (McFadden & McElhinny, 1990). The overall mean directions between the two lithologies are indistinguishable as they also pass a common mean test (Figure 8, McFadden & McElhinny, 1990). For the anorthosite sites that do not pass a common mean test with their diabase hosts, they nevertheless have coherent specimen-level directions that are close to their host diabase directions (Figure 8). We also plot the tilt-corrected mean pole of sites from both lithologies (diabase: 32. 5°N, 189.5°E, N = 15, A95 = 6.3, k = 37.4; anorthosite: 30.9°N, 190.8°E, N = 17, A95: 5.2, k = 48.5; Table. (Figure 8), consistent with the coeval magmatic activity between the Beaver River diabase and the PLV. If it is included in future Laurentia APWP compilations, it is this cooling unit mean pole paired with the estimated diabase emplacement age of 1091.7  E 0.2 Ma that should be used.

Thermal History Model
The consistency of the paleomagnetic directions between the anorthosite xenoliths and the host diabase indicate that the anorthosites were heated above the Curie temperature of low-titanium titanomagnetite (  E 580°C) within the Beaver River diabase. To determine whether this thermal history is consistent with the geometry of the units and to gain more insight into the emplacement history of the xenoliths, we developed a cooling model. In this model, the anorthosite xenoliths are considered to be solid spheres with an initial cool temperature embedded in a uniform sheet of diabase magma (Delaney, 1987;Unsworth & Duarte, 1979). The modeled thermal histories for various sizes of anorthosite xenoliths are shown in Figure 9. In one end member case, the initial temperature of the anorthosites is assumed to be 50°C. While this temperature is unrealistically low given that the anorthosites likely have a deep crustal source, thermal modeling shows that even a 100-m anorthosite xenolith with such low initial temperature would have been heated to the temperature of the tholeiitic magma (1150°C) within the sill. This temperature is well above the Curie temperature of magnetite. Anorthosite xenoliths with an assumed initial temperature of 500°C will equilibrate with the magma temperature on a similar, but slightly shorter, timescale. Therefore, the model predicts that the remanent magnetizations of the anorthosites will be reset during emplacement within the diabase sills, regardless of their initial temperatures. Model parameters set to match the xenolith AX16, from which a U-Pb date was developed in this study, leads to a model where the 27.5 m xenolith would have stayed at the magma temperature for about 100 years after sill emplacement (Figure 9). This duration estimate is a minimum as it does not consider heating associated with melt in the lower crust or during ascent prior to emplacement although this was likely rapid. The xenolith would have then cooled through the Curie temperature of magnetite (580°) after  E 1 kyr and acquired its magnetization as it cooled through magnetite blocking temperatures (down to  E 500°, Figure 7).

Origin and Age of the Anorthosite Xenoliths
There have been divergent interpretations regarding the age and magma source of the anorthosite xenoliths in the Beaver River diabase (Figure 1). Grout and Schwartz (1939) recognized the xenolithic nature of the anorthosites and suggested that the massive intrusion of the older anorthositic gabbro within the Duluth Complex may have supplied anorthosite fragments that were later entrained by the Beaver River diabase emplacement. Morrison et al. (1983), on the other hand, argued that the xenoliths were sourced from Paleoproterozoic or Archean lower crust that were liberated and contaminated by MCR magmas based on Sm and Nd isotopic data. They interpreted a Sm-Nd model age of 1.9 Ga from one of the xenoliths as providing a minimum crystallization age for the anorthosites though they acknowledged that these constraints are not definitive with respect to the age.
In contrast to this Archean to Paleoproterozoic model, Miller and Chandler (1997) favored a scenario where the anorthosite crystallized as part of MCR magmatism. They cited work by Kushiro (1980) who showed that the changing density contrast between labradoritic to bytownitic plagioclase and tholeiitic magma at different crustal pressures would promote flotation of plagioclase in deep (  E 20 km) crustal magma chambers and the creation of anorthosite cumulates in the lower crust. This mechanism of plagioclase flotation likely created massive anorthosite cumulates in the roof zones of subcrustal magma chambers during MCR magmatism. Miller and Weiblen (1990) speculated that plagioclase-phyric magmas tapped from these deep chambers fed shallow (  E 5 km) subvolcanic intrusions of the Duluth Complex, thereby creating the troctolitic anorthosites and gabbroic anorthosites of the Anorthositic Series. Miller and Chandler (1997) suggested that the nearly pure anorthosite xenoliths occurring in the younger and more hypabyssal diabase intrusions of the Beaver Bay Complex were harvested from these phase-segregated intrusions in the lower crust. They further argued that the isotopic data of Morrison et al. (1983) can be explained by anorthosite-forming MCR magmas having been contaminated by older crust rather than the anorthosites being older lower crust that was contaminated by MCR magmas.
Our new geochronology documents that the anorthosite xenoliths were liberated from depth and were emplaced within the shallow intrusions of the Beaver River diabase at 1091.7  E 0.2 Ma (95% CI). This timing of emplacement is constrained by the Beaver River diabase postdating the new 206 Pb/ 238 U zircon date of 1091.83  E 0.21 Ma for the AX16 xenolith and being older than the cross-cutting 1091.61  E 0.14 Ma Silver Bay intrusives.
The most straight-forward interpretation of the anorthosite xenolith 1091.83  E 0.21 Ma U-Pb zircon date is that it records crystallization of the anorthosite cumulates during Beaver Bay Complex magmatism just before the time of Beaver River diabase emplacement. The significant negative Eu anomaly in the zircons within the anorthosite constrains them to have crystallized from a magma that had experienced significant plagioclase extraction (Figure 6; Rubatto, 2002;Schaltegger et al., 1999). This result indicates that the zircons were comagmatic with their host anorthosite plagioclase. The Ti-in-zircon temperature estimates indicate that they crystallized from temperatures of  E 998 to 860°C (Supporting Information S1; Ferry & Watson, 2007). In addition, zircons that have lower Ti-in-zircon temperatures have lower Eu abundance, but enrichment of incompatible elements such as Hf and Th (Supporting Information S1). This systematic pattern of elemental concentration variation is consistent with the zircons crystallizing from residual melts on a cooling path that increased incorporation of incompatible trace elements and deepened the Eu anomaly with decreasing temperature and melt fraction. Scanning electron microscopy on two undated anorthosite xenoliths with plagioclase laths displaying interlocking textures shows zircon crystals with subhedral to anhedral shapes within the mineral assemblage that is interstitial to the plagioclase (Supporting Information S1). CL images show internal zoning in zircons which can be attributed to variations in REE, particularly Dy elemental concentrations, during zircon crystallization ( Figure 5; Remond et al., 1992). These data confirm that the zircons formed from residual melt within the interstitial spaces of the plagioclase cumulate and are inconsistent with a later metamorphic origin.
This scenario requires that there were large lower crustal magma chambers in which flotation of plagioclase resulted in cumulate formation during ca. 1,092 Ma Beaver Bay Complex magmatism and contrasts with the model of Miller and Chandler (1997) for an older origin of the anorthosite during the ca. 1,096 Ma Duluth Complex magmatism. Zircon U-Pb dates nearly always record crystallization age as the temperatures necessary for significant diffusive Pb loss exceed typical liquidus temperatures of zircon-bearing rocks. However, the anorthosites are a rather unique case given that the melting point of anhydrous plagioclase with an average composition of the Beaver River anorthosite (∼70% E anorthite, Doyle, 2016;Morrison et al., 1983) is quite high at  E 1400°C. Thermal modeling indicates that the anorthosite xenoliths would have equilibrated to the temperature of the olivine tholeiitic magma (  E 1100-1200°C) and remained at that temperature for more than 100 years in the diabase sill interior (Figure 9). While these temperatures would not have melted the plagioclase or zircon, they are high enough to consider the possibility of Pb diffusion out of zircon. Could diffusive resetting of the zircon in the anorthosite cumulates xenoliths allow their crystallization at ca. 1,096 Ma in the deep crust, but the closure of U-Pb zircon chronometer upon emplacement and cooling at ca. 1091.8 Ma?
The magnitude of Pb diffusion is dependent on the time spent at such a temperature. Using the diffusion parameters of Cherniak and Watson (2001), a sustained temperature of 1200°C for  E 10,000 years is required for diffusive loss of  E 90% of Pb from a  E 120  E m diameter zircon. In this case, zircons that crystallized at 1,096 Ma and then lost  E 90% of their Pb at 1091.6 Ma could give apparent U-Pb dates of 1091.8 Ma that are reproducible at the measurement resolution ( Figure 10). However, CL imagery reveals sharp boundaries between zones of differing CL response ( Figure 5) on the scale of  E 2  E m. Such CL zoning patterns are dominantly attributed to concentration variations in the REE Dy (Remond et al., 1992). A time-temperature history that results in 90% Pb diffusion out of a 120  E m diameter zircon would also cause Dy re-equilibration throughout a zircon, leaving no clear zonation ( Figure 10; Cherniak et al., 1997). Therefore, a scenario where the zircons first crystallized during Duluth Complex magmatism and subsequently lost more than 90% of Pb is exceedingly difficult to reconcile with the preservation of such thin, sharp zones. In fact, preservation of REE zoning in these zircons limits heating at the emplacement temperatures of the Beaver River diabase to a duration more consistent with our modeled duration of  E 100 years of heating prior to cooling to the temperatures that preserve such zonation (Figures 9 and 10). It is therefore most probable that the Beaver River diabase anorthosite xenoliths are entrained cumulate enclaves that formed at the time of Beaver Bay Complex magmatism.

A Comagmatic Relationship Between the Beaver River Diabase and the Greenstone Flow
Given the existence of many anorthosite xenoliths whose short-axis diameters often reach tens of meters and can be as wide as 180 m (Figure 1; Boerboom, 2004;, the Beaver River diabase magma conduits must have been at least this wide during magma ascent. It would be consistent with such wide conduits extending to hypabyssal depths for magma that flowed through these conduits to have vented to the surface. The high volume of the extrusive Greenstone Flow of the PLV leads to a potential match for this large feeder system. Doyle (2016) proposed a comagmatic link between the Beaver River diabase and the Greenstone Flow. Doyle (2016) discovered that both the intrusive Beaver River diabase and the Greenstone Flow have indistinguishable primary compositions that followed similar differentiation patterns. Doyle (2016) also highlighted the shared petrographic textures between the ophitic Beaver River diabase and the ophitic Greenstone Flow, which features the plagioclase laths clustering together and joining along their long crystallographic axes. The fosterite content of the olivines and enstatite content of the pyroxenes in the Beaver River diabase together with the Silver Bay intrusions, and the Greenstone Flow have overlapping compositions consistent with the same magma source (Figure 11). The composition of the plagioclase within the units further strengthens this interpretation. Although there are no known multi-crystalline anorthosite xenoliths in the Greenstone Flow, plagioclase megacrysts occur in the lava flow. Analyses of the anorthite content from plagioclase megacrysts show very similar m. Curves represent time-temperature conditions under which zircon will lose the indicated fraction of total Pb. Middle: Modeled zircon Pb loss scenarios with initial crystallization ages of 1091.8, 1,096, and 1,800 Ma with varying degrees of Pb loss at 1091.6 Ma compared to the actual U-Pb dates. Bottom: Conditions for preservation of Dy zoning in zircon. Curves represent time-temperature conditions under which different zoning thicknesses would be preserved in zircon. For conditions above the upper solid curves in each group, welldefined zoning will be lost at a given thickness. For conditions above the dashdot lines, zones will be partially lost but will retain initial composition in zone center. Pb diffusion and Dy zoning models follow Cherniak and Watson (2001) and Cherniak et al. (1997). Figure 11. Box plots of geochemical analyses of olivine, pyroxene, and plagioclase in the Beaver River diabase (BRD) and Greenstone Flow (GSF). The fosterite content in olivine crystals and the enstatite content in clinopyroxene crystals are very similar in the Beaver River diabase and the Greenstone Flow. The anorthite concentrations in the core, groundmass, and rim of the plagioclase megacrysts within the Beaver River diabase and the Greenstone Flow share very similar patterns and the distributions are nearly identical. The box encloses the middle 50% of the data ranges (i.e., the interquartile range), and the notch represents the median values. The whiskers extend to the 2.5th and 97.5th percentile values. An, anorthite; En, enstatite; Fo, fosterite. Data from Doyle (2016). values between the Beaver River diabase and the Greenstone Flow basalt (Figure 11; Doyle, 2016). In both units, the plagioclase cores are more enriched in anorthite than the rim and the groundmass. These data provide evidence that the core of the plagioclase megacrysts in the Greenstone Flow derived from a similar source with those in the Beaver River diabase and that the rims are later overgrowths. These mineralogical similarities are consistent with the interpretation that the Beaver River diabase and the Greenstone Flow have the same magma source.
The synchroneity between the Beaver River diabase and the Greenstone Flow inferred from comparable lithologies and geochemistry can be further evaluated using the paleomagnetic pole positions and radioisotopic dates from both units (Figures 4 and 8). The heat diffusion model of the cooling history of the anorthosite xenoliths within the diabase suggests that the time it takes to cool the diabase and anorthosite from low-titanium titanomagnetite Curie temperature (  E 580°C) to their blocking temperatures (  E 500°C) is on the time scale of a few thousand years (Figure 9). This time scale is close to the typical 4 10 E years which is considered to be sufficient for averaging out secular variations of the geomagnetic field. Figure 8 shows the site mean paleomagnetic pole positions from all diabase and anorthosite sites in this study against the previously synthesized Laurentia APWP developed using an Euler pole inversion to chronostratigraphically constrained volcanic poles in present-day coordinates (Swanson-Hysell et al., 2019). The site mean pole positions of the diabase and anorthosite overlap within uncertainty ellipses and the mean pole positions fall between the 1,095 and 1,090 Ma pole path positions (Figure 8), consistent with the geochronology results ( Figure 4). Further, the mean paleomagnetic pole position derived from the Greenstone Flow share a common mean with those of the Beaver River diabase and the anorthosite xenoliths, but these poles do not share a common mean with the mean pole derived from the PLV (Figure 8; Swanson-Hysell et al., 2019). This result suggests that the timescale over which the Beaver River diabase and the Greenstone Flow acquired their magnetization may be too short to fully average out secular variation. In this case, the overlapping pole positions between the Beaver River diabase and the Greenstone Flow strengthens their temporal correlation even more ( Figure 8).
The U-Pb dates are consistent with a comagmatic relationship as they reveal indistinguishable ages for the Beaver River diabase and the Greenstone Flow. The age of the Beaver River diabase is constrained to be between the 206 Pb/ 238 U dates of 1091.83  E 0.21 Ma and 1091.61  E 0.14 Ma (Figure 4) (Figure 4).
The PLV, including the Greenstone Flow, are interpreted to have erupted into the main central graben of the MCR during an interval of significant subsidence ( Figure 12; Cannon & Hinze, 1992;Miller & Chandler, 1997). In contrast to the thick accumulation in the PLV, the Beaver Bay Complex has an erosional (and slightly angular) unconformity atop it that is then covered by the younger Schroeder-Lutsen basalt (Figure 1;Miller et al., 2001). This relationship suggests that the Beaver River diabase was emplaced into a rift flank highland that experienced uplift during the active development of the central graben (Swanson-Hysell et al., 2019). Eruptions fed through the Beaver River diabase network would have emerged from the rift flank and flowed from the highland into main rift basin ( Figure 12).
The proposed intrusive-extrusive connection between the Beaver River diabase and the Greenstone Flow would imply that the Greenstone Flow extended for more than 250 km from northeastern Minnesota to the northern end of Isle Royale where the flow is  E 100 m thick and to the northeastern end of the Keweenaw Peninsula, where the flow is  E 400 m thick (Figure 1). With this length, the full volume of the Greenstone Flow reaches  E 6,000 3 km E (Doyle, 2016), rivaling the largest known lava flows on Earth (Figure 3). Such lengths were achieved for multiple high-volume flows within the Columbia River basalts (Reidel et al., 2013) and are modest compared to the Rajahmundry Trap lavas of the Deccan Traps which traveled  E 1,000 km (Self et al., 2008). One potential challenge for a flow having traveled from the Beaver River diabase to the PLV is that reconstruction of present-day rift basin isopachs from seismic data indicates that there is a deep bowl-shaped volcanics-filled basin offshore of Minnesota and the Beaver Bay Complex (Stewart et al., 2018). While this basin has a thicker accumulation of volcanics than surrounding regions, it is unclear whether it was a topographic barrier. If it was, it could have prevented lavas from present-day northern Minnesota from reaching the portion of the basin now exposed on the Keweenaw Peninsula. Therefore, it is also possible that the indistinguishable ages and similar geochemistry between the Beaver River diabase and the Greenstone Flow are the result of them having been derived from a contemporaneous deep magmatic source without being connected on the surface.

Conclusion
There was voluminous emplacement of magma into the shallow subsurface and eruption into the MCR basin ca. 1091.7 Ma at the end of the main stage of MCR volcanism. The anorthosite xenoliths within the Beaver River diabase and their U-Pb geochronology, whose interpretation is informed by REE patterns, indicate that there was a contemporaneous deep crustal magmatic system in which flotation of plagioclase formed anorthosite cumulates. The large dimension of the anorthosite xenoliths require that conduits feeding magma to the surface had widths that exceeded 150 m. These conduits would have delivered a high volume of magma into the rift basin. The high-precision U-Pb dates, together with paleomagnetic and geochemical data, are consistent with the hypothesis that the Beaver River diabase was the feeder system of the Greenstone Flow although they could have been disconnected at the surface and both be emblematic of this high-volume pulse of magmatism.