Geophysical Constraints to the Geological Evolution and Genesis of Rare Earth Element–Thorium–Uranium Mineralization in Pegmatites at Alces Lake, SK, Canada

Abstract

This investigation establishes an integrated method for rare earth elements (REE) exploration through a very promising and advanced exploration prospect in the Alces Lake area (SK, Canada) by assessing the integrated analysis of several multisource geophysical datasets. The resulting outcome provides important lithostructural information to the well-exposed, mineralized middle-to-lower crust at Alces Lake, comprising deep-seated poly-phase folds, ductile shear zones, and brittle faults. Geophysical–geological models of the Alces Lake property were constructed at different scales. The area of interest is located within the Beaverlodge Domain, about 28 km north of the Athabasca Basin’s northern margin. It contains some of the highest-grade rare earth elements (REE) in the world with the REE hosted predominantly in monazites within quartzo-feldspathic granitic to biotite–garnet–monazite–zircon-rich restite-bearing/cumulate mush melt pegmatites of anatectic origin (abyssal). Geophysical magnetic, gravity, and radiometric data were used together with Shuttle Radar Topography Mission (SRTM) images to facilitate the processing, modeling, and interpretation. Consequently, major structures were identified at different scales; however, the emphasis was given to studying those at the district/camp scale. The REE zones discovered to date occur within a large district-scale refolded synformal anticline. The eastern limb of this folded structure comprises a 30–40 km long, NW-trending shear zone/fault corridor with deep-seated structural crustal roots that may have served as the major pathway for ascending fluids/melts and facilitated the emplacement of mineralization. Thus, shear zones, faults, and folds in combination with lithological contacts/rheological contrasts appear to control residual/cumulate pegmatite emplacement and monazite deposition. Anomalies obtained from the airborne equivalent thorium survey data prove to be the most useful for REE pegmatite exploration. The results herein provide new interpretation and modeling perspectives leading to a better understanding of the distribution and lithostructural controls of REE on the property, and to new guidelines for future exploration programs at Alces Lake and elsewhere in northern Saskatchewan.

1. Introduction

The demand for the rare earth elements (REE) has increased substantially in recent years due to numerous new applications in a variety of manufacturing technologies, including high-tech, military, clean energy (e.g., wind turbines and permanent magnets), transportation, nanotechnology, medical equipment, etc. [1,2,3]. Suddenly more mining and exploration activities are required and more skilled earth scientists (in particular geophysicists with exploration capabilities) are needed to meet the growing demand [4]. While REE deposits appear to be abundant across the globe, economically viable deposits are comparatively rare. Since China, the major exporter of REE, restricted the export of this resource to other nations in 2010, many kinds of research and collaborative efforts have needed to be initiated to find alternative sources for the metals and develop new supply chains [3,5,6]. Granitic (abyssal) pegmatites, which are typically formed under upper amphibolite-to-granulite facies conditions and constitute products of partial melting [7], are one of the important sources of REE, which could help us achieve these goals. Thus, there is a resurgence of interest in exploration for rare earth elements in pegmatites across Canada, in particular the Precambrian Shield areas, to ensure the availability of REE to the different technology industries (e.g., [3,8]). Normand [8] has documented the extensive pegmatite occurrences/fields in northern Saskatchewan that are mineralized with REE, including the areas surrounding the prolific uranium-bearing Athabasca Basin of Mesoproterozoic age. We acknowledge that the mineralized pegmatites at Alces Lake are unique, enigmatic, and controversial in their origin.

Rare earth pegmatites are important sources of many critical elements including rare metals, uranium (U), thorium (Th), tantalum (Ta), niobium (Nb), tin (Sn), lithium (Li), beryllium (Be), and boron (B) [9]. The pegmatite fields/occurrences of northern Saskatchewan are comparable on the basis of the nature of their enrichment to other typical deposits of this type located in Russia (Aldan, Khibina Massif), Sweden (Ytterby), Canada (Tanco), among others [10,11,12,13,14,15,16,17]. Even though many of the rare metals used today are found within rare metal granitic pegmatites, debate continues about the origin of these rocks; especially those located within the lower-to-middle crust (i.e., abyssal type) [7,11]. Understanding the mineral systems of the rare metals is an integral part of making new discoveries and thus reducing supply deficits. Successful exploration programs depend on the ability to access and apply fundamental knowledge around mineral systems, deposit type descriptions, and ore genesis [18].

During the exploration programs, many methods are used in order to discover new REE deposits. This study illustrates the most effective steps to follow to define the mineralization bodies and potential drilling targets (i.e., provides an exploration template). Due to the presence of abundant thorium–uranium-bearing minerals in abyssal mafic-rich pegmatites, radioactivity is a common characteristic of this type of REE mineralization, with the radiation being primarily a result of the ubiquitous presence of thorium in the REE minerals (e.g., monazite and/or allanite). So, apart from using handheld radiation detectors, such as Geiger counters and scintillation counters (scintillometers) and portable XRF instruments, geophysical methods (including both ground and airborne geophysical techniques) can also be used to detect and map potential prospects/ore bodies. Various geophysical data are collected, processed, and analyzed at different scales, thus an integrated geophysical prospecting approach can contribute to new discoveries on a regional/province- to district/camp- to deposit/prospect scale. Even though many of these geophysical methods (except for airborne and ground radiometrics) cannot provide direct detection of REE mineralization, they can use complementary processing tools, coupled with 2D/3D modeling, to provide valuable information/insights about some of the geological/lithological controls. This in turn can guide us to identify structural pathways and/or traps of the given mineralization under investigation [19,20,21].

Herein, we document the geological–geophysical modeling research of the low- to high-grade REE mineralization at Alces Lake, which is hosted by REE-rich pegmatite bodies. In order to characterize the pegmatite field of the Beaverlodge Domain for its REE potential, Appia and the current researchers have carried out geochemical analyses of the REE-bearing and associated crystalline rocks of Archean/Paleoproterozoic age, as well as geophysical surveys to locate and map the mineralized pegmatite bodies. However, currently ongoing exploration efforts did not effectively utilize the geophysical characterization of the potential mineralization zones in the Alces Lake area by advanced data processing, interpretation, and integration of multi-source geophysics. This is considered necessary to delineate geologic structures and associated geochemical traps that could provide sites to the REE mineral deposits/prospects within the area.

Presently, the prominent mapping tool is gamma-ray spectrometry surveying for REE mineralization at Alces Lake. Gamma-ray spectrometry datasets essentially reflect geochemical signatures of potassium (K), uranium (eU), and thorium (eTh), which are related to the mineralogy and geochemistry of rocks underlying the region, as well as any overlying soil/regolith. Prior to this study, successful implementation of all potential gravity data were lacking (i.e., the different derivatives, as demonstrated in this investigation).

Thus, here we identify and document the detailed structures of the surface/subsurface using various tools to analyze and interpret the aeromagnetic data, gravity, aeroradiometric, and Shuttle Radar Topography Mission (SRTM) imagery data. Potential field measurements of the Earth’s crust and upper mantle provide important information on the subsurface structures of crystalline basement rocks [22]. The results of potential field mapping provide information on alkaline igneous complexes and middle crustal anatectic pegmatite provinces that are often related to pegmatite REE mineralization (e.g., [22,23,24,25,26,27,28] among others). However, the application of integrated geophysical methods to map the subsurface geology in terms of structural fabrics/trends, presence of hydrothermal alteration zones, identification of granite intrusion (source rock), and elevated uranium (eU) and thorium (eTh), which are related to REE mineralization in the Alces Lake area, is still evolving both in terms of more datasets, more advanced processing, and more integration to understand the minerals system at different scales (i.e., regional/province to district/camp to deposit/prospect).

The objective of the present study was to map and assess REE mineralization in the pegmatite-rich Alces Lake area by using ground geophysics, ground magnetics, aeromagnetics, and airborne gamma-ray spectrometry data with DEM/LIDAR survey data. In this manuscript, we first present the regional geological setting of the Beaverlodge area, followed by the deposit-scale geology and a brief outcrop/petrological description of rocks exposed within the study area. In the Data and Methodology section we present the datasets and methods that have been used during the processing and interpretation stages. The Interpretive Results section follows. The manuscript ends with the Discussion and Conclusions sections, where we review the advantages and disadvantages of the proposed methods and compare obtained results with other research. By integrating the new geophysical datasets and drill hole data with existing geological and structural data of regional-to-deposit scale (i.e., [29,30,31] among others), we produced detailed geophysical maps and models of the subsurface related to the main showings of the Alces Lake property. The study also applies the geophysical characteristics of the main WRCB–Ivan–Dylan–Dante area of known low- to high-grade REE abyssal granitic pegmatites in the area to assess the mineralization potential of other underexplored/unexplored areas of the larger property area. Thus, ultimately, our contribution provides a better understanding of the geological evolution and structural controls of low- to high-grade REE mineralization (i.e., the minerals system—see ref. [18]) within the entire Alces Lake area.

2. Regional Geology

Alces Lake is located approximately 34 km east of Uranium City and about 28 km north of the Athabasca Basin margin. It is sited within the Beaverlodge Domain near the junction of the Beaverlodge, Train, Zemlak, and Ena domains in the Rae Subprovince (Figure 1 and Figure 2). Beaverlodge Domain is characterized overall by lower magnetic signatures than the Train Lake Domain and Zemlak Domain, and high gravity signatures (Figure 2).

The Beaverlodge Domain consists of polydeformed ortho- and paragneisses. It is separated from (1) the Train Lake Domain by the Oldman-Bulyea shear zone (OBsz) and the Alces Lake shear zone (ALsz)/St. Louis fault (SLf), (2) the Zemlak Domain by the Black Bay fault (BBf), and (3) the Tantato Domain by the Grease River shear zone (GRsz) (Figure 1, Figure 2 and Figure 3). The rocks identified throughout the Beaverlodge Domain include ca. 2.6–3.1 Ga high-grade Archean orthogneisses and migmatitic paragneisses, ca. 2.33 to <2.17 Ga Murmac Bay Group, ca. 1.83 to 1.818 Ga Martin Group, and ca. 1.75 to 1.50 Ga Athabasca Group [32]. Two structural levels have been identified by regional mappers of the study area, and separated by what was interpreted as a folded unconformity surface [32,33].

Rocks of the lower structural level comprise 3.4–2.6 Ga high-grade orthogneisses and paragneisses, which are widespread in other parts of the Rae craton [34]. These rocks record the effects of Neoarchean (ca. 2.57 Ga) and Paleoproterozoic (ca. 2.5–2.3 Ga, 1.94–1.90 Ga, 1.90–1.85 Ga, 1.84–1.80 Ga) tectonism of the Arrowsmith, Taltson, Snowbird, and Trans-Hudson orogenies. At this structural level, the following lithostratigraphic units have been identified: 3.06–2.99 Ga granites, 2.68–2.60 Ga granites, and 2.33–2.29 Ga granites, found in association with migmatitic paragneisses [33]. There are also ca. 1.94–1.90 Ga structures, which are associated with E-SE-striking transposition foliations and tight-to-isoclinal folds that were refolded into NW-trending open folds [33]. They have been overprinted by NE-striking fabrics, folds, and dextral mylonitic shear zones associated with dextral transpression along the Grease River shear zone at ca. 1.90–1.85 Ga [33,35,36,37].

Rocks of the upper structural level are supracrustal rocks of the 2.33 to <2.17 Ga Murmac Bay Group [33], a strongly metamorphosed sedimentary succession, consisting of basal conglomerate and quartzite to psammite that is locally intercalated with dolostone, iron formation, and mafic volcanic flows, all of which are overlain by psammopelite to pelite gneisses [32]. The rocks of the Murmac Bay Group have been affected by greenschist- to granulite-grade metamorphism, with the peak metamorphic conditions being attained during 1.94–1.90 Ga tectonism, which are the result of burial and accompanying lateral compression, thrusting, and strike-slip shearing during the Taltson orogeny [33].

Two other lithological packages/groups occur regionally outside the investigated area; the Martin Group in the Uranium City area and the Athabasca Group south of the study area (see

Table 1. Template for mineral exploration showing the main lithostratigraphic units within the immediate study area. Wavy line represents the unconformity.

). In the western Beaverlodge Domain, the Martin Group lies unconformably on top of the Murmac Bay Group and its underlying basement rocks, and comprises continental red beds as well as mafic rocks. The rocks of the Martin Group have not been metamorphosed, but have undergone open folding and faulting [38]. To the south and west, the Athabasca Group package of rocks unconformably overlies the Murmac Bay Group rocks and locally the Martin Group rocks. This package comprises flat-lying, unmetamorphosed sedimentary rocks that are composed mainly of quartz-rich sandstone of MesoProterozoic age [39,40].

Within the Beaverlodge Domain, the mineralization is closely associated spatially with the contact between Archean rocks and Paleoproterozoic supracrustal packages composed of quartzite, amphibolite, psammite, psammopelite, and pelite rocks (Figure 1 and Table 1) of the Murmac Bay Group [33,41,42].

3. Alces Lake Geology and Mineralization

Alces Lake is ranked by grade for containing the highest-grade REE occurrences in Canada [41] and the second highest in the world after the Gakara Rare Earth Elements (REEs) deposit [43]. Currently, the Alces Lake tonnage is unknown, as there is no preliminary NI-43-101 resource estimate. The REE in most of the high-grade zones are hosted mainly in monazites (locally the monazite content reaches up to 85 modal %,) within quartzo–feldspathic granitic to biotite–garnet–monazite–zircon-rich restite-bearing/cumulate mush partial melt pegmatites. The REE mineralogy of the low- to medium-grade zones (low-grade mineralization is 0.1 to 1.0 wt.% TREO (total rare earth oxides), medium-grade mineralization is 1.0 to 4.0 wt.% TREO, and high-grade mineralization > 4.0 wt.% TREO) is currently under investigation by Appia Rare Earths & Uranium Corporation of Toronto, Canada (formerly Appia Energy Corp. and abbreviated to Appia hereafter in this paper). The REE-bearing pegmatites are associated spatially with regional- to district/camp-scale fold structures, ductile shear zones, and brittle faults, with the latter two structures cutting both paragneiss and orthogneiss [8,42,44]. The Alces Lake property is located within the Beaverlodge Domain (northern Saskatchewan) about 28 km north of the Athabasca Basin margin [29,33] (Figure 1 and Figure 2). The mineralized pegmatites show a similar structural setting to the Fraser Lakes and Kulyk Lake pegmatites [45,46] that are localized along or near major deep-seated crustal shear zones/faults, immediately east of the eastern Athabasca Basin [47,48].

The Alces Lake REE mineralized area is a part of a regional refolded fold structure, within a synformal anticline. The mineralization is located on the eastern limb, close to the hinge of a south-plunging, truncated open fold [8,42] (Figure 3, Figure 4 and Figure 5). It is enriched with critical REE (neodymium (Nd), praseodymium (Pr), dysprosium (Dy), terbium (Tb) used for permanent magnets. At a 4.0 wt.% total rare earth oxide (TREO) cutoff for the high-grade zones, the Alces Lake average grade is 16.65 wt.% TREO [49].

Lithological units on the Alces Lake property include:

• Archean granitic gneiss;
• Paleoproterozoic metasedimentary gneiss (pelitic-psammopelitic [+/−graphite], quartzite, amphibolite, pyroxenite, diatexite), and feldspathic gneiss of the Murmac Bay Group;
• Paleoproterozoic syn-to late-anatectic pegmatites;
• Paleoproterozoic late-orogenic-to-metasomatic biotite schist, pegmatite augen, and monazite accumulations (the REE mineralized system) [49].

Table 1 provides more detailed descriptions of the main stratigraphic and lithological units within the study area.

The REE mineralized system is composed of Proterozoic late-orogenic to metasomatic massive braided biotite schist, quartzofeldspathic pegmatite augen, and monazite accumulations. In outcrop, the mineralization is unique and enigmatic. Most of the REE in the high-grade mineralized zones (e.g., the Wilson and Ivan zones) are hosted primarily within monazites, and to a lesser degree in allanite, xenotime, and apatite, within both the biotite schists (glimmerites) and pegmatite augens [49,52]. The monazite mineralization occurs in different habits from isolated grains within 1 to 3 cm thin lenses, as isolated massive clusters up to meters thick, and as pervasive massive clusters of monazite in augenitic (boudinaged) masses along ductile to ductile/brittle to brittle structures. The authors opine that the mineralized pegmatites have been emplaced within/near the Archean/Paleoproterozoic transition zone under middle crustal P-T conditions and form polyphase anatectic pods/boudins/zones along/near this transition [53]. Additionally, emplacement of the pegmatites is interpreted to be deep-seated, structurally controlled along ductile to ductile-brittle to brittle shear zones/faults associated spatially with regional polyphase megafold structures [53,54].

The classification of the mineralized pegmatites and associated rocks follows that of Cerny and Ercit [7], after Cerny et al. [55]. Here, below in Figure 6, we present some of the field characteristics and lithological associations of these pegmatites, which are integral to understanding their geophysical characteristics and using different exploration geophysics tools in their discovery from surface to depth; i.e., the focus of this paper.

4. Data and Methodology

Lithostructural and modeling interpretations of the Alces Lake area, using various multisource geophysical datasets, were combined with SRTM imagery data.

High-resolution airborne magnetic and radiometric surveys, and drone magnetics and ground gravity, were conducted for the study area in 2016 and 2019, respectively. The geophysical datasets and accompanying data processing reports were provided by Appia Energy Corp. Mostly these types of airborne surveys are used for the exploration purposes in order to detect certain types of mineralization, such as iron oxide–copper–gold deposits, massive sulfides, and heavy mineral sands, kimberlites, skarns, porphyritic intrusions, etc. These geophysical surveys can also be an effective means to map faulting, alterations, and basement geology beneath the sediment cover (such is the case, for example, of unconformity-related uranium mineralization). Furthermore, SRTM imagery can be very effective while carrying out the structural and lithological interpretations of the studied area, contributing a lot to the understanding of the mineralized zones as well.

Table 2. Types and sources of data used for geophysical analysis.

Data Type	Files	Source
Maps and Figures	Geological map of the Alces Lake property	Appia Rare Earths & Uranium Corp.
	Detailed figures of the Alces Lake zones
Reports	Alces Lake Field Assessment Report, 2013
	Alces Lake Report, 2013
	Report on Airborne VTEM, Magnetic and Radiometric Survey, 2016
Geophysical data	Drone Magnetic Survey Data, 2019
	Ground Gravity Survey Data, 2019
	Magnetic and Radiometric Survey Data, 2016
	Regional data (TMI, Local phase of the residual total magnetic field; Bouger gravity)	[56]
Topography	DEM	Appia Rare Earths & Uranium Corp.
	Topography data
	SRTM data (Shuttle Radar Topography Mission)	[57]
Geological data	Bedrock geology (ArcGIS shapefiles)	[31,58]
	Faults (ArcGIS shapefiles)
Worm Data	Saskatchewan Geological Survey Open File 2010-48	[59]

and Figure 7 provide the types and sources of data and the workflow, respectively, which were used for geophysical analysis in this research investigation. Additionally,

Table 3. The main types of input data with examples and rationale in their use for REE exploration purposes.

Scale	Type/Source of Data	Geophysical Methodology Used	Data Obtained as a Result	Rationale	Direct Identification of REE Pegmatitic Bodies
Province/ Regional	Magnetic data	TMI -> RTP -> 1VD/tilt-angle -> lineament features	Delineation of major geological structures, such as faults, shear zones, some lithological boundaries.	Vertical and tilt-anglederivatives of magnetic data reveal responses from near-surface sources, which help identify the spatial locations of major geological structures.	No
District/Camp to Deposit/Prospect	Magnetic data	TMI -> RTP -> 1VD/tilt-angle -> lineament features	Delineation of local scale geological structures, such as faults, shear zones, some lithological boundaries, Archean–Proterozoic transition zone.	Shallow, higher-order, local faults can serve as structural traps and also provide a pathway for mineralizing fluids. Such features are represented as linear features on geophysical data. Vertical derivative and tilt-angle of magnetic data enhance responses from near-surface sources; thus, the identification of these geological structures is facilitated.	No
	SRTM, Topographic data, Google Earth images	Lineament analysis	Delineation of local scale geological structures, such as faults, shear zones, some lithological boundaries.	Topographic data help mapping the surface morphology and relief of the ground. Identification of major faults, shear zones, and some lithological boundaries.	No
	Equivalent Th data	Elevated values of eTh	Delineation of REE-Th pegmatites and/ortheir geochemicalhalos.	Radiometric data reflect the surface geology, and due to the correspondence between the REE and Th within the Alces Lake deposit (observed during the mineralogical studies), we can see the corresponding radiometric eTh anomalies over the mineralized areas.	Yes
Deposit/Prospect	Magnetic data (drone-borne)	TMI -> RTP -> 1VD/tilt-angle -> lineament features	Delineation of deposit-scale geological structures, such as faults, shear zones, some lithological boundaries.	Drone-borne and ground-based magnetic data allows modeling of magnetized buried geological bodies with higher precision.	In part
	Gravity data	Bouguer gravity -> Residual -> 1VD/tilt-angle/Apparent density -> lineament features	Identification of hidden ore bodies using the density contrast among rock units. Mafic (monazite-biotite-zircon +/− garnet–rich) pegmatite bodies have higher gravity values.	Allows the delineation of the outlines of lithological units and monazite-biotite-zircon +/− garnet–rich pegmatite bodies (provided a large enough density contrast exists between the mafic pegmatitic bodies and their host rocks).	In part, sometimes yes

gives the main types of input data with examples and rationale in their use for REE exploration purposes.

The geophysical data (observed gravity and magnetic fields) show the combined effect of all the sources from multiple depths, so it is important to differentiate between them. Thus, different filtering techniques were implemented in order to separate the anomalies of different wavelengths from one another. The shallow sources anomalies are often represented by short wavelengths (or high frequencies) and the deep source anomalies correspond to long wavelengths (or low frequencies) [60].

4.1. Geophysical Methods

Methodologies of this investigation for all phases of the geophysical data, their advantages, and disadvantages followed the currently prevailing applications of geophysics in mineral exploration [61]. The effectiveness, capabilities/limitations, and costs of all methods are compared (see Figure 1.2 of Dentith and Mudge [61]). Resolution and adaptability of the most advanced data management, data processing, and interpretation techniques are documented in detail. Within the present investigation, an example of the successful advanced application of the filtering technique is given below in the district camp-scale magnetics section, where the results of high-pass filtering clearly map the images of a shallow complex fold structure, with wavelengths of around 250 m, indicating that the source is not more than ~200 m depth and that it extends to ~1000 m depth.

Regional magnetic data used in the current research investigation were obtained from the Geological Survey of Canada (GSC) geophysics portal in Ottawa, Canada. The GSC regional surveys were flown at different parameters (dependent on vintage), some at 300 m altitude and others at 90 m altitude, with flight line spacing varying from 3200 m to 400 m resolution, respectively. The authors acknowledge that these GSC spacing sets are not sufficient for identification of the Alces Lake mineralization package of rocks at shallow depths; thus, tighter 200–100 m spacing would give much better resolution and interpretation ability. Additionally, regional magnetic and gravity multi-scale edges (worms) from the Saskatchewan Geological Survey (SGS) were also used during the analysis, and represent a combination of horizontal derivatives and upward continuations of geophysical data. These worm data were evaluated to delineate and interpret the different geological variables at corresponding depths.

To address resolution and interpretation needs, district/camp- to deposit/prospect-scale geophysical gravity and magnetic surveys were carried out by Appia for its area of interest. Magnetic data are a cost-effective way to acquire geological information, as well as being an effective structural mapping tool [61,62].

The potential field datasets were provided by Appia to the first author. This field information was collected by Geotech Ltd. Of Aurora (Aurora, ON, Canada) and MWH Geo-Surveys Ltd. Of Vernon (Vernon, BC, Canada), respectively. The airborne magnetic survey was conducted from 28 to 30 May 2016, with a total of 147 line-kilometers of data being acquired at a mean altitude of 74 m above the ground [63] (Figure 8 and Figure 9). A more detailed outcrop-scale drone survey was conducted on 8 June 2019. These two surveys were done at different scales: the Geotech survey was heliborne at 100 m spaced lines, whereas the MWH drone magnetics survey was acquired at 25 m spaced lines (at ca. 30 m above the ground) and the MWH ground gravity was acquired at 10 m spaced lines. Gravity station spacing was 4 m × 10 m. Figure 8 shows the boundary outlines of the geophysical surveys at Alces Lake.

4.2. Aeromagnetic Data

Magnetic anomaly images (Figure 8) show amplitude variation of approximately 275 nT for the district/camp-scale area and 130 nT for the deposit/prospect-scale area due to the variation of susceptibility values of different lithological units. The general regional patterns/trends of the magnetic field are NNW-SSE, E-W, and NE-SW. The magnetic lithologies of the study area include a magnetite-bearing mixture of Fe-rich paragneiss including silicate facies iron formation, amphibolite, pyroxene-biotite-garnet gneiss, and garnet diatexite, whereas non-magnetic lithologies are represented by quartzitic to psammopelitic to pelitic gneisses and migmatites with ubiquitous intrusions of S-type, ilmenite-bearing leucosomes, pegmatites, and leucogranites. Geophysical and geological data reveal that the investigated strongly peraluminous pegmatite bodies are essentially nonmagnetic at/near the surface.

The observed magnetic field responses represent the combined response from the sources of variable depths, where shallower sources often produce the shorter-wavelength anomalies and deeper ones produce the longer-wavelength anomalies. So, in order to obtain the magnetic data that best represented the studied geological structure/body, several data enhancement techniques were used to selectively improve the signal of interest [61,62], and were applied to the ‘raw records’ data, in conjunction with field observations. Such techniques allow geological structures to be observed and mapped in considerably more detail than utilization only of the ‘raw records’ data.

The IGRF was removed and reduction-to-pole operator [61,62] was applied by mathematically transforming the observed magnetic field to the equivalent field that would have been observed at the magnetic north pole (assuming purely induced magnetism; i.e., no remanent magnetism) (Figure 10).

Detailed structures can be revealed based on the study of the first vertical derivative (1VD = ∂B/∂z, where ‘∂B/∂z’ represents the partial derivative of the field with respect to depth (z), Figure 11d and Figure 12d). It is typically used in order to accentuate the potential field expression of shallow sources [61,62]. The 1VD is applied to the data by taking the numerical derivative of the data with respect to the vertical (depth) direction. The result of this operation is a transformed set of data that highlights subsurface features that have a strong vertical gradient.

Tilt-angle derivative (TDR) was used to highlight the shallow features [61,62] by incorporating horizontal as well as vertical gradient information (TDR = arctan(VD/THG), where THG is the total horizontal derivative (gradient), and VD is the vertical derivative) (Figure 11 and Figure 13). These two enhancement techniques show the inferred faults as well throughout the study area. The 1VD image shows that magnetic linear features along the NW direction are dominant in the eastern part of the area and along the NE direction in the northern area (Figure 11d). The inferred faults around the more magnetic body in the southwestern part are also more visible (Figure 11d).

High-pass and low-pass filtering are techniques used in signal processing to remove unwanted frequencies from a signal [61,62]. These techniques are utilized to see deeper (low-pass) or shallower (high-pass) features of the subsurface.

Total horizontal gradient (THG) determines both shallow and deeper structures [61,62]. The THG is a two-dimensional derivative that calculates the horizontal gradient of the data by computing the gradient in both the x and y directions and then obtaining a vector sum (THG = sqrt(Hx × Hx + Hy × Hy), where Hx and Hy are cross-line and in-line horizontal gradients). It reflects the distribution and lithological contacts of the rocks of the area, showing the contrast between the more magnetic garnetiferous diatexite or S-type granite and less magnetic psammopelitic to pelitic gneiss, migmatite, and diatexites (Figure 11b).

4.3. Gravity Data

A high-resolution ground gravity survey was conducted for the study area by MWH Geo-Surveys Ltd. Beginning 4 June and concluding 18 June 2019 [64]. The data were collected at 1495 stations. The gravity survey was carried out at 4 m × 10 m station spacing (Figure 14); considered to be leading edge in terms of the typical spacing for gravity surveys. However, parameters such as outcrop and downhole density measurements of the different lithological units within the grid area are also considered to be as important for constraining the gravity models [65]. Gravity data were first processed using the standard methods (e.g., ref. [66]) followed by the tilt-angle, 1VD, etc. (Figure 15 and Figure 16), and residual gravity (Figure 14b), which ultimately allowed us to carry out the apparent density calculations (Figure 15e and Figure 16c).

Apparent density mapping transforms Bouguer gravity anomalies into density distribution [61] (L(r) = r/2πG(1−e^−tr), where L(r)—apparent density, G—gravitational constant, t—thickness of the earth model). For this transformation, first the 1000 km Butterworth high-pass filter was applied to cut off the long-wavelength signal. Then, the filtered dataset were used for apparent density calculation. We set the average crustal density to 2600 kg/m³ and defined the thickness of model layer to be 200 m. This apparent density map shows the lithological units well, as there is a strong correlation between the rock types and apparent density. The areas with the highest density are those where there are mafic pegmatite bodies containing high-density monazite mineralization, zircon, biotite, garnet, etc. (Figure 16d).

All of the filters and enhancements that we applied give us a much better understanding of the distribution of rock types.

4.4. Radiometric Data

Radiometric data were acquired by Geotech in 2016 together with the aeromagnetic data [63]. Radiometric surveys have proved themselves to be very useful during surface prospecting and geological mapping, representing the upper 30–100 cm of the subsurface. Radiometric data give approximations for potassium, uranium, and thorium (K, U, and Th) concentrations. It is noted here that Th and U concentrations are normally expressed as equivalent Th and equivalent U [27].

Airborne radiometric surveys, combined with magnetic and gravity surveys, are routinely used in various exploration programs because they provide an opportunity to better show, understand, and interpret the subsurface and its rocks, in particular for the rocks that produce significant radiometric anomalies, such as different types of gneisses, K-rich feldspathic rocks, and U- and Th-rich rocks. The lithological contacts can be best delineated/discerned from the K and eTh maps. Presence of the monazite-rich pegmatitic rocks shows high Th contents (Figure 17). The K map (Figure 17b) shows similar results.

4.5. SRTM Images Analysis

SRTM images are very useful for structural and lithological discrimination of the study area and are used extensively to identify the tectonic/structural regime, which contributes as well to increasing our knowledge on the particular types of mineralization [67,68]. Apart from identifying the individual lineament features, they can also be useful to conduct a lineament density analysis to evaluate the structural complexity of the area. They can assist to identify major intersection zones of possible fault dilation with associated fluid and heat flow. These are the places where structural geochemical traps and/or pathways for mineralizing fluids can be found. Thus, gradient S, N, W, E, NE, NW, Laplacian filters, etc., were applied and various vector layers were created that represented faults, shear zones, and fold axes for further use and interpretation.

5. Interpretive Results

The geological framework of the Alces Lake area is recorded in the regional-to-district to prospect-scale magnetic, gravity, and radiometric data and on outcrop-scale magnetic and gravity maps. While regional geophysical datasets reflect the orientation and position of province-scale anomalies and major geological units and structures, interpretive analysis was made possible with the use of available deposit/prospect and outcrop-scale datasets. This allows a much closer look at small-scale geological structures in our area of interest. With the help of the regional magnetic and gravity multi-scale edges (i.e., worms—cf., [69]), we were able to recreate the spatial location of major geological structures; notably the Black Bay Fault and the Alces Lake Shear Zone/St. Louis Fault (Figure 1). The major synclinal (synformal) form of the district-scale area was also traced from the worms’ data up to several kilometers in depth. It is recognized that due to the nature of the worms, their upward continuations become less reliable and less interpretable with increasing depth [69]. Regional magnetic and gravity maps showed a strong correlation between the magnetic and gravity characteristics and the major geological structures; notably, faults and tectonic domains. Thus, it was possible to trace the contours of the Beaverlodge Domain, which has low-magnetic and high-gravity signatures, and the shear zones/faults that separate it from the Tantato Domain (by the Grease River shear zone (GRsz) and more magnetic Train Lake (by the Oldman-Bulyea shear zone (OBsz) and the Alces Lake shear zone (ALsz)/St. Louis fault (SLf)) and Zemlak domains (by the Black Bay fault (BBf)).

The magnetic data were not successful in showing the direct expression of the REE (Figure 13). However, commonly the edges of magnetic anomalies can help to locate geological contacts such as lithological contacts, fault zones, and shear zones. Additionally, it is possible to detect whenever there exists a significant contrast between the magnetization properties of two geologic bodies [28]. Thus, aeromagnetic data in combination with SRTM imaging and regional magnetic data revealed the location of most of the major geological structures observed/interpreted at the regional-to-district scale via the manual lineament feature analyses method (Figure 3, Figure 4 and Figure 5). Several stages of tectonic evolution can be mapped this way including different generations of folds, ductile shear zones, and ductile–brittle to brittle faults. Using this methodology, the main synformal anticline, the Alces Lake shear zone—St. Louis Fault and late cross-structures were delineated, and as well importantly from an exploration perspective, a major shear zone lying immediately east of the main synclinal form of the area was mapped (Figure 3, Figure 4 and Figure 5). Younger brittle faults that were mapped at this scale were overall EW-trending. With the use of SRTM images, dominant lineament features were identified, occurring along NW-SE, NE-SW, and EW directions. The brittle EW faults are interpreted to be the latest post-tectonic features.

The present integrated gravity and magnetic study delineated the subsurface structures as well as their association with the mineralized system. By studying the district/camp- and deposit/prospect-scale magnetic data, the authors were able to distinguish and differentiate between several rock units (however, it has to be mentioned that the deposit/prospect-scale magnetic data contributed less to the final interpretation and the district/camp-scale data were mostly used). In general, the observed magnetic highs represent a mixture of paragneiss, pyroxene-biotite-garnet gneiss, and garnet diatexite and magnetic lows by quartzitic-to-psammopelitic to pelitic gneiss and migmatite (Figure 12 and Figure 13). Similarly, gravity highs within the outcrop-scale data are associated with the garnet-rich rocks, amphibolites, glimmerites, and high-grade monazite mineralization, and lows with various quartz-rich rocks, including Murmac Bay metasediments and non-mineralized pegmatites (Figure 16).

Radioelement surveys and especially radioelement ratios can be a highly effective means of directing mineral exploration efforts [25,62]. For the studied area, it is suggested that the eTh maps are the most useful for REE exploration due to their significant number of stronger responses, many of which are associated with mineralized bodies. Thus, the equivalent thorium (eTh) image proved to be extremely helpful for targeting the monazite-rich pegmatites (and biotite schists) over the area of study, as very strong correlations were found between the radiometric maps and locations of monazite-hosting rocks (Figure 17). The majority of the strong eTh anomalies lie (generally but not exclusively) within the lows of more magnetic units—paragneiss, pyroxene-biotite-garnet gneiss, and garnet diatexite. The highest eTh concentrations seem to be localized in proximity to the hinge zones of major fold structures. However, it is important to recognize that relying solely on one survey image or geophysical method may not necessarily provide success. Nevertheless, incorporating radiometric maps in the exploration process for pegmatite-hosted rare metals is considered highly beneficial and is suggested to be used for the other areas of interest/potential on the Alces Lake property.

6. Discussion

Geophysical data can provide much insight into an individual deposit. Utilizing this data make it possible to characterize the geological structures that may have influenced the emplacement of the REE mineralized bodies and create detailed geophysical maps and models of the subsurface of the studied area. With ever-increasing demand for a great variety of minerals, in remote complex areas, concerted efforts are being dedicated to develop airborne, mainly aerial vehicle (drones-based) instrumentation [70]. Additionally, software development is being extended in several areas: (1) modification and integration of existing programs to process and permit the combined interpretation of large volume and diversified datasets (e.g., cloud-based PCA and K-means cluster analysis), (2) incorporation of deep learning (AI) principals to new approaches in the analysis of different but interrelated data volumes [71], and (3) development of simultaneous real-time modeling procedures for complex data volumes. There are numerous examples of successful applications of the gravity, magnetic, and radiometric survey data for research and exploration purposes for rare metals [23,24,25,72]. For example, a gravity survey conducted at the deposit scale at Tanco deposit [25,73] was successful at direct delineation of the pegmatite-hosted rare metal mineralization due to the large density contrast between the host and mineralized rocks. Magnetic data can be used for both direct identification of metal deposits (e.g., identification of carbonatite deposits, which often correspond to highly magnetic anomalies of circular shapes [25,74] or delineation of major structures, which can serve as the pathway and/or trap for their emplacement ([53,54,75], this research).

The methods and techniques used in the current study can be applied as well to the rest of the Alces Lake property and to similar prospects within the Beaverlodge Domain to increase our overall understanding of this REE-rich domain and to provide a template for future exploration. For example, airborne magnetic data have the ability to identify the main structural trends within the area and radiometric data can be used to locally target the REE mineralization on the surface directly. Hence, more detailed ground gravity surveys can help outline the shape and depth of individual bodies or packages of rock (Table 3). However, it has to be stated that relying solely on one survey image or geophysical method may not necessarily provide successful results and thus careful consideration of other existing geological information is essential for correct interpretations. It is also important to keep in mind the limitations of each of the methods. For example, some of the main methods (like radiometrics) show a major decreasing resolution with depth (i.e., different for each type of survey) and/or non-uniqueness of interpretation. It is also worth mentioning that the Alces Lake mineralized area is unique and the exploration workflow suggested in this paper would best work for analogous mafic-rich pegmatites.

The filtered magnetic anomalies identified for the Alces Lake deposit/prospect suggest that the location of the mineralized zones is structurally and lithologically controlled. Such interpreted maps also show that mineralization is localized in proximity to the linear magnetic anomalies that represent the axially planar fold nose zones or the sheared limbs of major fold structures. This supports an exploration model that is premised on the fact that REE mineralization took place along structures during and after their formation (i.e., syn- to late- to post-deformation). Integration of SRTM and geophysical information facilitated identification of lithological and structural units within the study area. As a result of the above-documented modeling exercises and the application of various enhancement/filtering techniques, the configuration and extension of geological layers were confirmed at Alces Lake, and new potential geological bodies were identified based on their respective geophysical characteristics. By combining the information from different sources, we have interpreted the geophysical data, produced the resultant constrained profiles and models, and generated an interpreted geological map from geophysics (Figure 18), which can better show/explain the location of mineralized units with respect to geological units and major structures. Additionally, major lithostructural controls on the location of mineralization have been identified at various scales from lineament analysis of magnetic/gravity/radiometric images superimposed on geological maps; from district to local to deposit scale.

The mineralized pegmatitic rocks and biotite schists (glimmerites) at Alces Lake are interpreted by Poliakovska et al. [53,76] to have formed by partial melting and subsequent fractional crystallization during the Taltson-Thelon Orogen event. The geochemical data show that the pegmatites are predominately peraluminous with a lesser volume being metaluminous, which suggests that the pegmatite-forming melt probably had a source similar to the migmatitic metasedimentary gneiss units within the exhumed lower to middle crust of the Alces Lake area, with possibly some contribution from Archean orthogneisses and mafic-rich mantle rocks as recorded by the dense (gravity high) Beaverlodge Domain (Figure 2). These mineralized rocks are probably similar to those of the Fraser Lakes area [45], which are classified as abyssal middle crust anatectic pegmatites [7,11,55]. This interpretation means that such pegmatites have an anatectic origin and thus are not directly related to a parental granitic pluton. They are expected to be overall mineralogically simple and to lack any significant zonation [45]. However, this is not the case at Alces Lake, whereby the mineralogy becomes complex as a result of hot restite–rich melt emplacement, country rock interaction/reaction, and crystallization (see detailed description below). Overall, their bulk composition depends on the composition of the metasedimentary and associated rocks that underwent anataxis at upper amphibolite to granulite facies conditions [77].

The OBsz hosts the east-northeastern boundary of the Beaverlodge Domain (Figure 1). It is interpreted to be a shear zone/reactivated thrust fault [78], which wraps around to the NW and then westwards to become the Alces Lake Shear Zone/St. Louis Fault, and to the SE and then abruptly SW to become the Grease River Shear Zone. Thus, the base of the Beaverlodge Domain is interpreted to be a décollement surface (verging ENE) that dips moderately and then shallowly to the west-southwest. This is supported by geological map observations and interpretation in the eastern-most Beaverlodge Domain, whereby when one approaches the OBSZ, outcrops of dense mafic gneiss/granulite (lowermost Beaverlodge Domain) become prevalent. The latest dextral sense of movement on the Grease River Shear Zone is dated at circa 1855 Ma [36]. Hence, this surface is interpreted by the authors as being the crustal level where mafic magma underplating (i.e., ponding of dense hot mafic magmas) took place, and then followed in sequence by crustal heating, partial melting, magma mixing, and generation of mafic-rich residual melt/cumulate rocks. Various large-scale structures (weakened shear/fracture zones) facilitated crustal melt transfer. A number of mafic dykes (now amphibolites, interpreted as metamorphosed mafic rocks of mantle origin) of different generations (and ages) are mapped within major structures throughout the Beaverlodge Domain and along its boundaries, in particular along/adjacent to its northern boundary, the Alces Lake Shear Zone [51].

The mineralized zones are interpreted to represent a shear/fracture zones along which the melts/fluids were introduced. Overall, the rocks of the area have undergone dynamic metamorphism, which is expressed by ductile protomylonitic to mylonitic fabrics, including augen structures, mafic boudins within the amphibolite units, stretching and thinning of leucocratic material, and the presence of leucocratic boudins in the melanocratic gneiss matrix. The mylonitic fabrics are ubiquitous throughout the paragneisses. Based on field observations, the authors suggest that the presence of glimmerites (biotitite)–biotite-rich rock is the result of the interaction/crystallization of a ‘hot’ residual/cumulate magma with Fe/Mg/S-rich metasediments that resulted from the upward migration of partial melting products from the source area [79]. According to the model described by Silva et al. [80], glimmerites can form as a result of melt flow ascending using the pathways in high-strain zones, leading to melt-rock interaction (i.e., high-T metasomatic mineral reactions). Such deeply seated structures/shear zones are considered major pathways for melts originating via anatectic processes. During reactivation pulses, melt phases can migrate through these rheological weakened parts of the crust and penetrate along the network of fractures, shear zones, and lithological contacts to facilitate/enhance the processes of host rock assimilation and metasomatic mixing, thus leading to the formation/crystallization of glimmerites. Similarly, such a process is proposed here for many of the mineralized zones in the Alces Lake area. Biotite schists and biotite-rich rocks are found extensively in the vicinity of the shear zones along with the restite-rich, cumulate mush, monazite-rich pegmatites, and associated quartzites, paragneisses, and amphibolites.

In summary, the mineralization at Alces Lake is interpreted to be plausible restitic and cumulate melt material formed from ubiquitous disequilibrium melting at depth because of interpreted anomalous heat from mafic magma underplating. As a result, an extremely restite-rich melt phase migrated to the presently observed crustal level at Alces Lake (i.e., old middle crust at upper amphibolite/lower granulite facies) (Figure 19).

7. Conclusions

Multifaceted geophysical techniques have to be crucial components of the exploration strategies for REE deposits. Magnetic data can help identify the main structural trends of the area, while radiometric methods can be used to target the REE hosts. The ground gravity survey data (and 100 m spaced airborne gravity surveys) can consequently be employed to conduct more precise exploration when a target has been identified. Therefore, having access to a wide range of airborne data (e.g., magnetic, radiometric, electromagnetic, and gravity) greatly improves the likelihood of discovering new economic prospects.

An attempt to analyze the structural setting of the Alces Lake area through the geophysical, SRTM, and radiometric datasets was highly successful. This was achieved by utilizing several enhancements and filtering techniques, including 1VD, TDR, THG, HP, LP, tilt-angle, and apparent density.

Detailed geophysical survey data (i.e., both airborne and ground) from the Alces Lake property area provided insight into known and prospective areas of REE mineralization. The magnetic field data show the distribution of the magnetic properties (signature) within the rock units and highlight the major structures and lithological units within the area of study, which in turn increases our understanding of the complexity of the geologic and tectonic history of this area. Our template presented here for rare earth element (REE) exploration provides the methodology for integration, interpretation, and application of multisource geophysical datasets within a highly complex geological environment and tectonic framework hosting REE mineralization.

This investigation supports the geological model in which the mineralized pegmatites were emplaced within/near the Archean/Paleoproterozoic transition zone with REE being hosted within polyphase anatectic pods/boudins/zones (residual/restite-cumulate mush-rich crustal melts) along/near this transition. Lithological and structural contacts were demarcated using a combination of SRTM images and airborne geophysical data. These observations also indicate the very high potential utilization of gravity surveys at proper station/line spacing in the search for REE mineralization. The principal structures identified for the study area occur along NW-SE, NE-SW, and EW directions, where the brittle EW faults have been interpreted to be the latest post-tectonic features. These results are in agreement with and add to the structural trends deduced from other geological sources. Large-scale structures and their flexures played a big role by controlling the emplacement of mineralization at different structural levels. These structures probably served as a major transcrustal pathway for ascending crustal melts/fluids and facilitated the emplacement of mineralization. Thus, large-scale folds, shear zones, and faults, in combination with lithological contacts/rheological contrasts, are interpreted to control restite-bearing/cumulate mush-rich pegmatite emplacement and extensive enigmatic monazite deposition at Alces Lake.

REE-hosting rocks are characterized by high-Th anomalies, dense minerals, and para- to non-magnetic mineral assemblages. Thus, it was possible to successfully recognize the corresponding radiometric high-eTh anomalies over the mineralized areas by radiometric data correlation with the surface geology, and associated high correspondence between the REE and Th concentrations within the Alces Lake mineralized zones (observed in mineralogical studies).

This comprehensive investigation also reveals that a systematic planning of an orderly application of geophysical methods from beginning to end is necessary to minimize cost and effective development of a prospect/deposit. In addition, future studies at Alces Lake would benefit from a better characterization of the polyphase structures throughout the study area, which would improve our understanding of the overall lithostructural controls related to the REE mineralization. Finally, detailed petrological, petrophysical, and geochronological studies are needed to improve our knowledge of the genesis of the REE mineralization and the timing(s) of emplacement.