Controls on, Expressions of and Discovery Potential for Gold Mineralization in the Central-Eastern Yilgarn Craton, Western Australia: New Insights from an Integrated Targeting Study

1. Introduction

The central-eastern Yilgarn Craton of Western Australia (Figure 1), subject to this study, hosts a gold endowment of >40.0 Moz Au [1] that is mainly contained within 11 orogenic gold deposits, or deposit clusters, of Archean age. Given this significant endowment, the presence of 78 operating gold mines and the potential for new discoveries, interest in the area has been considerable, both in terms of gold exploration [2] and gold-related research activities [3,4,5,6,7,8,9,10,11,12,13].

Despite the economic importance of the region, published targeting studies for orogenic gold deposits in all, or part, of the central-eastern Yilgarn Craton are scarce, limited to the works of [14,15,16,17,18]. In 2000, Groves et al. [14] established a comprehensive conceptual targeting model that served to inform a knowledge-driven vectorial fuzzy logic approach [19] to mineral potential modeling (MPM), covering an area of ~17,000 km² in the southern Kalgoorlie Terrane (as such, spatial overlap with our study area is minimal). The model was based on only two types of input data: The known gold occurrences and a lithostructural (solid) geology map. The authors interpreted the MPM results in terms of gold mineralization being controlled by only a handful of critical factors such as proximity to NW-SE-striking granitoid-greenstone contacts, crustal-scale faults, in particular NW-SE-striking fault segments, NW-SE-striking lithological contacts, proximity to anticlinal hinges, and chemically-reactive rocks.

In contrast to [14], a 2010 study by Czarnota et al. [15], which covered both the Kalgoorlie and Kurnalpi terranes (~130,000 km²), benefited from a decade of progress in computing power, significant advances in the capabilities of geographic information systems (GIS) as well as MPM, a vastly improved understanding the geology of the Yilgarn Craton and its orogenic gold deposits, and the emergence of the mineral systems concept [20,21,22], which introduced a shift in targeting by moving away from descriptive, deposit-centric constructs to a broader, scale-integrated approach focused on the processes of mineral deposit formation and their geodynamic context. In contrast to the earlier study by [14], Czarnota et al. [15] developed numerous geological, geochemical and geophysical predictor maps, believed to represent mappable features of the targeted gold systems. However, rather than using mathematical algorithms (i.e., MPM) to combine these predictors, the authors employed what was described as a “simple additive GIS treatment.” Unfortunately, this pseudo-MPM approach prevented the use of spatial statistics to assess the validity and performance of the input predictor maps and the model as a whole. Instead, [15] defined an arbitrary buffer zone around areas of high gold prospectivity, which they claimed captured 75% of the known gold endowment in 5% of the modelled area.

The most recently published MPM study by Witt et al. [16] developed gold prospectivity maps for both the Eastern Goldfields Superterrane (EGST) (~235,000 km²) as well as the Yilgarn Craton as a whole (~680,000 km²). The authors quantified spatial relationships between known gold endowment and a comprehensive set of regional geological, structural, geochemical and geophysical criteria, from the best of which they built prospectivity maps utilizing a knowledge-driven fuzzy logic approach [23,24]. The best predictors of gold potential at the regional scale, as highlighted by the authors’ spatial statistical analyses, were as follows: Proximity to mafic granite intrusions, areas of high fault density, regional fault bends, and greenstone belts.

This MPM-assisted targeting study, undertaken over an area of ~27,150 km², differs from the earlier work in three key ways: (i) It is geophysically driven, using enhanced gravity and magnetics to map structural architecture, including cryptic basement structures; (ii) it uses enhanced radiometrics and remotely sensed image data to map prospective alteration systems; (iii) it creates a new, internally consistent lithostructural framework that integrates those datasets into a single solid geology base for MPM; and (iv) it applies and cross-validates five complementary MPM techniques, spanning knowledge-driven, data-driven and continuous methods. Together these steps allowed us to reduce the search space by an order of magnitude and identify underexplored gold prospective corridors largely blind to traditional outcrop-biased methods.

2. Materials and Methods

2.1. Study Background

The study described in this paper took place in late 2021. It was commissioned by Australian mining company St Barbara Limited, which also determined the layout of the study area. At the time, the central-eastern Yilgarn Craton was home to St Barbara Limited’s Leonora operations, consisting of the Gwalia underground mine and 1.4 Mtpa processing plant, as well as nearby gold development and exploration opportunities. In 2023, the Leonora operations were sold to fellow Australian gold miner Genesis Minerals Limited [25].

2.2. Methodology

The approach taken in this study entailed six principal steps:

• A review of the orogenic gold deposit model with a focus on deposits of this type within the study area and broader Yilgarn Craton.
• The generation of a mineral systems-type targeting model [20,21,26,27,28], emphasizing the critical processes required to form orogenic gold deposits and their mappable expressions.
• The filtering of geophysical and remote sensing data and extraction of structural information and hydrothermal alteration signatures from these data [29,30].
• A solid geology (lithostructural) interpretation based on all available geological, structural and geophysical information.
• Mineral potential modeling (MPM), adopting a multi-technique approach that entailed the use of continuous as well as data- and knowledge-driven mathematical techniques and, thereby, facilitated the cross-validation and comparison of the resulting prospectivity maps [31].
• The identification of gold exploration targets.

2.3. Supporting Data

The geoscience and exploration data used in this study had two principal sources (Table 1): Open-access repositories maintained by the Geological Survey of Western Australia (GSWA) and proprietary company data.

2.4. Processing and Interpreation of Geophysical and Remote Sensing Data

2.4.1. Geophysical Data and Processing

Geophysical data were downloaded in grid format (as per Table 1) and windowed to the study area, allowing for adequate buffer zones. Minimal processing (preparation of the data for filtering) was required as the data were available as merged products with the appropriate data reduction complete, providing blanket coverage of the study area at various spatial resolutions (i.e., cell sizes). Access to objective, ‘blanket coverage’ geophysical data of adequate spatial resolution is essential for exploration targeting in areas where prospective bedrock is masked by significant cover like in the central-eastern Yilgarn Craton. In addition, geophysical data are ideally suited to mapping expressions of mineralising processes and, therefore, are a fundamental component of a mineral systems-type targeting approach.

Three types of geophysical data were sourced and filtered for inclusion in this targeting study [34]:

• Gravity data, which respond to variations in rock density and reflect density contrasts in the subsurface. As such, gravity surveys are effectively a ‘depth penetrative’ method.
• Magnetic data, which respond to variations in rock magnetism, are mostly controlled by magnetic susceptibility. Like gravity data, magnetic data reflect subsurface contrasts.
• Radiometric data, which respond to surface variations in naturally occurring gamma radiation, mostly originating from radioactive isotopes of potassium (K), thorium (Th) and uranium (U).

In the EGST, gravity gradients show a clear association to faults that host gold deposits [35]. The extraction of gravity gradients from gridded gravity data has been carried out in various ways since the method was first published in the late 1970s [36,37]. Here we adapted the phase congruency approach of [38], using the differential upward continuation operator of [39] instead of log Gabor filters, resulting in a multi-scale feature detection filter perfectly tailored for potential field (gravity and magnetic) data. Topological features [40] extractable from gridded potential field data at multiple scales are:

• Edges (i.e., maximum gradient curves), which represent traces of physical property boundaries across gravity or magnetic data grids, arising from structural breaks such as faults or shear zones, or lithological contacts.
• Ridges (i.e., surface maximum/peak curves), which correspond to ridgeline curves across gravity and magnetic data grids. Gravity ridges, for example, reflect peak gravity responses. In granite-greenstone terrain, gravity ridges often coincide with the thickest parts of dense, mafic rock units and, thus, work well in terms of mapping structurally thickened mafic-ultramafic rocks in the central parts of inverted greenstone basins.
• Valleys (i.e., curves of surface minimum), which correspond to valley line curves across gravity or magnetic data grids. Magnetic valleys, for example, often coincide with faults or fault corridors along which primary magnetite was partly, or mostly, destroyed due to fluid-rock interaction.

Radiometric data, on the other hand, work well in terms of detecting and mapping outcropping hydrothermal alteration systems. For example, studies by de Quadros et al. [41], Herbert et al. [42] and Shebl et al. [43] served to demonstrate the merit of this technique for identifying areas of gold-related white mica (sericite) alteration, using various computations involving normalized radiometric data and their ratios.

2.4.2. Remote Sensing Data and Processing

ALOS digital elevation data [44] were downloaded in grid format (as per Table 1) and treated in a similar fashion as described for the geophysical data above. Topographic elevation data reflect the shape of the Earth’s surface across the study area, where surface highs often represent units more resistive to weathering. As such, even in ancient peneplained cratonic environments such as the Yilgarn Craton, topographic data can be useful for mapping areas of out- or subcropping bedrock, structure in exposed bedrock (and structural trends therein), or past and present mining activities. Surface topography can also reflect geological processes, which, in turn, can inform exploration targeting, especially in well-exposed neotectonic environments.

Multi-spectral Sentinel-2 data were downloaded as multiband grids (as per Table 1) and windowed to the study area. The area has several salt lakes that can generate false positive results in some of the processing. These lakes were manually traced and removed from the data prior to processing. Indexes were calculated to highlight ferric-iron-bearing phases like goethite (Band4/Band2) and clay (Band11/Band12). Additionally, the iron feature depth (IFD) was calculated to highlight possible weathering sulfides in the regolith [45]. The three indexes (goethite, clay, and IFD) were combined to create a CMY ternary image. The image was converted to CEI L*a*b* color space and the Euclidean distance to black was calculated [46]. Black areas in the CMY image are where all three indexes have high values. These are the most likely areas to have exposed clay alteration assemblages and weathering sulfides. The distance to black was converted to an inverse distance so that high values in the distance grid were areas where all three indexes have relatively high values. The inverse distance grid was thresholded and polygonised.

2.4.3. Solid Geology Interpretation

The workflow adopted in our interpretation of the geophysical data and generation of a solid geology (lithostructural) map of the study area was broadly similar to that described by [47]. Namely, (i) constructing form lines, (ii) identifying magnetic rock units, (iii) characterizing geophysical domains, (iv) integrating the geophysical interpretation with other relevant information (e.g., geological maps), (v) identifying and defining structural elements, (vi) mapping and interpretating lithology, and (vii) generating an integrated lithostructural framework.

2.5. Mineral Systems Concept

The targeting model developed in this study was generated according to the mineral systems approach [20,21,22,26,48] and in a manner described in more detail by [27,30,49,50]. In short, the mineral systems concept views mineral deposits as localized expressions of a series of geological processes operating at different scales of time and space:

• Source processes that extract the essential mineral deposit components (i.e., melts and/or fluids, metals, and ligands) from their crustal or mantle sources;
• Transport processes that drive the transfer of the essential components from source to trap regions via melts and/or fluids;
• Trap processes that focus melt and/or fluid flow into physically and/or chemically responsive, deposit-scale sites;
• Deposition processes that drive the efficient extraction of metals from melts and/or fluids passing through the traps;
• Preservation processes that act to preserve the accumulated metals through time.

Where one or more of these critical processes fail to operate, a mineral deposit cannot form or be preserved. The probabilistic philosophy at the core of the mineral systems concept lends itself well to MPM [27,49,50] given the latter runs on mathematical algorithms. Another key strength of the mineral system concept is that it provides a robust yet flexible framework for formulating a holistic, process-based targeting model and for observing, mapping, and/or querying in the available geoscience data the expressions of the critical processes of mineral deposit formation [22,51].

2.6. Fry Analysis

Fry analysis [52] presents a simple geometrical approach to visualizing patterns in spatial point data. The method works by recording the distance and bearing from each point in a given dataset to every other point, either manually or computationally. Hence, for n points in a dataset, Fry analysis generates n² ‒ n spatial relationships. These point translations can then be visualized in autocorrelation diagrams, known as Fry plots, and rose diagrams, serving to enhance subtle trends in point patterns and reveal directions of maximum continuity of mineral deposit alignment [53,54,55].

2.7. Mineral Potential Modeling (MPM)

MPM, first introduced in the late 1980s alongside the emergence of geographic information systems (GIS) [56,57], has since matured into a robust, cost- and time-efficient tool for exploration targeting that also possesses a strong capability for big data analytics. More specifically, MPM can concurrently handle, integrate, process, and model the typically different and often very substantial geological, geophysical, geochemical, remote sensing, and drilling data developed and utilized in mineral exploration. It is also a tool that is well-adapted for effective screening and target generation within large search areas, be it a particular mineral district or belt or an entire country or continent. In general, MPM encompasses the following principal stages [21,51,58,59,60,61,62]:

• Genetic model stage: Recognition the critical geological processes responsible for the formation of the targeted deposit type to establish a conceptual deposit model.
• Targeting model stage: Translation of the genetic model into a targeting model in which the critical ore-forming processes are expressed through mappable targeting criteria (also referred to as spatial proxies, predictor maps, predictors, or targeting elements).
• Mathematical model stage: Assignment of weights to integrate the diverse predictor maps utilizing mathematical algorithms.
• Target identification and prioritized stage: Mapping and prioritization of the most prospective areas.

The mathematical modeling can be broadly subdivided into knowledge-driven, data-driven, hybrid data- and knowledge-driven, and continuous techniques that utilize logistic functions [24,58,63]. The selection of a suitable weighting approach is strongly governed by data availability, specifically the number of known mineral occurrences in support of the targeting model, commonly referred to as prospect locations. For instance, ‘brownfields’ regions are typically well explored and, thus, more data-rich, especially in the vicinity of known mineral deposits. For such regions, the weighting is commonly data-driven. Un- or underexplored ‘frontier’, ‘grassroots’, or ‘greenfields’ regions, on the other hand, are typically data poor and may contain only few, if any, mineral occurrences of the target type. In such search spaces, a knowledge-driven approach is often required, with MPM reliant on expert opinion [24,58]. In contrast, continuous weighting approaches require neither prospect locations nor expert opinion, with continuous spatial evidence offering superior predictive capability compared to discretized evidence [63].

In this study, we adopted a multi-technique approach to MPM [31] comprised of continuous (fuzzy gamma: [64]; geometric average: [65]; improved index overlay: [66]), data-driven (random forest (RF): [67]) and knowledge-driven (best worst method–sum additive weighting (BWM-SAW): [68,69,70]) modeling techniques.

3. Geology

3.1. Eastern Goldfields Superterrane (EGST), Yilgarn Craton

3.1.1. Overview

The study area is situated in the central-eastern Yilgarn Craton, a large, poorly-exposed Archean (~3730 to 2660 Ma) piece of crust comprised of extensive granite-greenstone and high-grade metamorphic gneiss terrains [71] (Figure 1). More precisely, the study area straddles the boundary between the Kalgoorlie and Kurnalpi terranes, which form part of the ~2960 to 2660 Ma EGST [7,12,72,73,74], one of the most important and best endowed Late Archean metallogenic provinces in the world [8]. Each of the terranes that make up the EGST (i.e., from west to east, the Kalgoorlie, Kurnalpi, Burtville and Yamarna terranes) is a distinct, structurally-bound entity characterized by major differences in stratigraphic ages and whole rock and isotopic geochemistry of its volcanic and intrusive rocks. [72,75]. The individual terranes can be further subdivided into lithostructural domains defined by internally consistent stratigraphy and structural history and bounded by interconnected fault systems interpreted largely from geophysical data [75].

3.1.2. Geology and Structure

The greenstone belts and major fault systems within the EGST typically strike N-S to NNW-SSE. Deep seismic reflection data portray the principal fault systems, which can be traced in gravity and magnetic data for up to 700 km along strike, as E-dipping, crustal-scale structures intersecting the Mohorovičić discontinuity. The supracrustal volcanic and sedimentary successions that define the greenstone belts are preserved in up to seven kilometer-thick synformal basins with greenstone stratigraphies generally dipping and facing away from granite domes flanking or coring the belts [7]. As summarized by [10,76], greenstone sequences across the EGST show a broadly coherent stratigraphy, typically comprising (from the bottom to the top) of:

• A basal volcanic package dominated by ~2720 to 2690 Ma komatiite and basalt;
• A sedimentary package dominated by ~2690 to 2670 Ma deep marine siliciclastic and volcaniclastic sedimentary rocks; and
• A ~2665 to 2655 Ma sedimentary package of fluvial and deep marine sedimentary rocks, referred to as ‘late basin’ fill, that rest unconformably on the older greenstone successions [77,78].

Older, ∼2960 Ma to 2770 Ma greenstone successions, mainly comprised of intermediate to felsic and mafic-ultramafic volcanic rocks, have been identified in the Burtville Terrane. These rocks have lithological and temporal affinities with the Youanmi Terrane in the western Yilgarn Craton [7].

Granitoids are abundant and widespread in the EGST, and can be classified into five groups [7,8,12,33]:

• High-Ca granites: Constitute the most abundant group of granites (∼60% of all EGST granitoids) and occur both within and external to greenstone belts. They are typically granite, granodiorite or trondhjemite in composition. Emplacement ages fall into three clusters at ∼2800 Ma, 2740 to 2650 Ma and, most commonly, 2685 to 2655 Ma.
• High field strength element (HFSE)-enriched granites: Make up 5 to 10% of all Yilgarn Craton granites. They are typically granite or, to a lesser degree, granodiorite in composition and occur mostly internal or marginal to greenstone belts. Emplacement ages fall into two clusters at >2720 to 2665 Ma and, most commonly, 2700 to 2680 Ma.
• Mafic granites: Make up 5 to 10% of all Yilgarn Craton granites. They are typically granite, granodiorite, tonalite, trondhjemite or diorite in composition and occur mostly internal or marginal to greenstone belts. A common spatial association with gold mineralization has been noted. Their emplacement ages fall into a relatively narrow time bracket from >2720 to 2650 Ma, although younger examples likely exist.
• Low-Ca granites: Are the second most abundant granite group by volume (∼20%) and mostly external to greenstone belts. Compositionally, the granites of this group are typically granitic or granodioritic. Their emplacement ages fall into a relatively narrow time bracket from 2655 to 2630 Ma. Low-Ca granites are considered products of partial melting of high-Ca granite source rocks.
• Syenitic granites: Are the least abundant granite group (<5%) and typically internal or marginal to greenstone belts. Granites of this group are syenite or quartz syenite in composition. A spatial association with gold mineralization has been noted. Emplacement ages cluster at 2650 Ma and 2655 to 2645 Ma.

3.1.3. Deformation History

According to [7], a combination of poor rock exposure and strain partitioning make it difficult to correlate structural information across the EGST. As such, its structural event history remains poorly constrained and subject to ongoing debate. The structural framework study by [7] addressed these difficulties using an integrated approach that involved the collection of >10,000 new structural observations and structural, stratigraphic, magmatic, metallogenic and geophysical data correlation and integration. Based on the results of this work, the authors proposed a six-stage (D₁ to D₆) deformation history (Figure 2):

• D₁: ENE-WSW-directed extension at ~2720 to 2670 Ma marked by rifting and greenstone deposition.
• D₂: ENE-WSW-directed shortening at ~2670 to 2665 Ma marked by the termination of volcanic activity and generation of NNW-SSE-trending upright folds and N-S-to NE-SW-striking dextral strike-slip and reverse faults.
• D₃: NE-SW-directed extension and extensional doming (i.e., core complex formation) at ~2665 to 2655 Ma marked by deep crustal exhumation and the formation of late basins with the late basin sequences constituting the first record of the deposition of granite detritus in the EGST.
• D₄: ENE-WSW-directed shortening at ~2655 to 2650 Ma (D_4a) marked by the tightening of earlier formed fold structures, WSW-directed thrusting along NNW-SSW-striking faults as well as the generation of NNW-SSE-trending upright folds and reverse faults; and WNW-ESE-directed shortening at (D_4b) marked by reactivation of and sinistral transpression along earlier formed NNW-SSE-striking faults and generation of ENE-WSW-striking thrust faults, which recorded transport to the NW and SE.
• D₅: NE-SW-directed shortening at ~2650 to 2635 Ma marked by dextral strike-slip movement along N-S- to NNE-SSW-striking faults and thrusting along NNW-SSE- to NW-SE-striking faults.
• D₆: Low-strain vertical shortening and horizontal extension at <2630 Ma marked by the development of crenulations.

3.1.4. Metamorphic History

The EGST recorded metamorphic conditions ranging from sub-greenschist to granulite facies developed over a ~130 m.y. time period from ~2750 to 2620 Ma. A broadly conformable Neoarchean metamorphic history is indicated by regional metamorphic isograds extending across terrane and domain boundaries [11,32]. The latter identified five regional metamorphic events across the EGST, which have been temporally and spatially integrated with the above structural event framework of [7] (Figure 2).

• M_a: Early-formed, low-P/high-T upper-amphibolite to granulite facies assemblages of the M_a metamorphic event are rare. They are restricted to magmatic arc-related, ~2730 to 2810 Ma greenstone sequences of the western Burtville Terrane and HFSE granites and ~2675 to 2715 Ma greenstones of similar affinity in the Gindalbie Domain of the Kurnalpi Terrane.
• M₁: This metamorphic event, which occurred at ~2750 to 2700 Ma, produced high-P/moderate-T assemblages preserved in relatively narrow, upper-amphibolite grade domains along major, crustal-scale fault zones. The structural setting, depth of burial and rapid exhumation of these rocks are consistent with partial burial of buoyant magmatic arcs, or terranes, during magmatic arc accretion events in subduction-like environments.
• M₂: This ~2680 to 2670 Ma, low-P/moderate-T metamorphic event produced a contact metamorphic pattern linked to the emplacement of large volumes of high-Ca granite melt into the upper crust, which were generated by partial melting of a down-going slab subducted beneath the Kalgoorlie and Kurnalpi terranes. The M₂ event coincided in time with the cessation of volcanism and D₂ crustal shortening.
• M₃: The M₃ metamorphic event is believed to have been caused by lithospheric extension at ~2665 to 2650 Ma, linked to the cessation of subduction with sag on the previously subducted plate causing slab rollback and extension of the overriding plate. The M₃ event coincided in time with the D₃, that is metamorphic core complex formation.
• M₄: Low-P/high-T metamorphism at ~2650 to 2610 Ma was likely triggered by delamination of the lower crust, resulting in mantle upwelling and the arrival of a thermal anomaly in the upper crust associated with widespread and voluminous low-Ca granite magmatism.

3.1.5. Geodynamic Implications

The available stratigraphic, structural, geochemical, and isotopic data support an interpretation of the Kalgoorlie, Kurnalpi and Yamarna terranes of the EGST as failed intracontinental rift basins, developed on older ‘proto-Yilgarn’ crust. The rift basins were subsequently closed, uplifted, deformed, and metamorphosed during the ~2670 to 2630 Ma Kalgoorlie Orogeny [74,80].

3.2. Kalgoorlie-Kurnalpi Rift

The ~2720 to 2600 Ma Kalgoorlie-Kurnalpi Rift (Figure 1), recognized by [12] based on stratigraphic, geochronological and isotopic data correlations across the Kalgoorlie and Kurnalpi terranes of the EGST, is believed to have formed in response to a major mantle input event that exploited weakened crust along the eastern margin of the Yilgarn proto-craton, comprised of the South-West, Narryer and Youanmi terranes. According to [12], the genesis and evolution of the intracontinental Kalgoorlie-Kurnalpi Rift is consistent with spatially and temporally overlapping plume-related magmatism in the Kalgoorlie Terrane and west-directed subduction to the east of the Burtville Terrane. This geodynamic model, however, implies that plate tectonic processes similar to modern-style plate tectonics operated as far back as the Archean eon, a contentious hypothesis that is subject to ongoing debate [81,82,83].

3.3. Kalgoorlie Terrane

The Kalgoorlie Terrane (Figure 1 and Figure 3), the westernmost lithostructural element of the EGST, is bound by the Ida-Waroonga fault system to the W and Ockerburry fault system to the E. The terrane is divided into 10 domains, five of which (i.e., Boorara, Jundee, Moilers, Ora Banda and Wiluna) occur within the study area. The terrane is mainly composed of 2710 to 2660 Ma greenstone and overlying, coarse clastic (‘late basin’) sequences dated at 2658 to 2655 Ma. Older (>2730 Ma) lithostratigraphic sequences appear to be restricted to the Norseman, Boorara, and Wiluna domains. Due to widespread cover, clearly defined lithostratigraphic sequences have only been established for the southern Kalgoorlie Terrane. Here, greenstones can be divided into the ~2710 to 2690 Ma Kambalda Sequence (tholeiitic and komatiitic mafic-ultramafic rocks) and the overlying, ~2690 to 2660 Ga Kalgoorlie Sequence (felsic volcaniclastic and epiclastic rocks with subordinate lavas). The late basin sequences are mainly composed of poorly sorted, polymict conglomerates that grade upward into sand- and siltstone. Kalgoorlie Terrane greenstones form highly deformed, curvilinear belts that are separated by variably deformed and metamorphosed granitoid rocks, intruded principally between 2760 and 2620 Ma [7,72,84,85,86].

3.4. Kurnalpi Terrane

The Kurnalpi Terrane (Figure 1 and Figure 3), which abuts the Kalgoorlie Terrane to the east, is bound by the Ockerburry fault system to the W and . The terrane is comprised of seven domains, five of which (i.e., Edjudina, Gindalbie, Laverton, Menangina and Murrin) occur within the study area. Lithostratigraphically, the Kurnalpi Terrane greenstones can be divided into four sequences: (i) The ~2800 Ma Laverton sequence (mafic-ultramafic volcanic rocks and banded iron formation with minor clastic sedimentary rocks intruded by felsic porphyries); (ii) the ~2715 to 2705 Ma Kurnalpi Sequence (mafic volcanic and intrusive, intermediate calc-alkaline igneous and feldspathic sedimentary rocks), which has no obvious stratigraphic equivalent in the Kalgoorlie Terrane; (iii) the ~2695 to 2675 Ma Minerie Sequence (mafic-ultramafic volcanic rocks), which are correlated to the Kambalda Sequence in the Kalgoorlie Terrane; and (iv) an unnamed, ~2690 to 2680 Ma sequence of bimodal rhyolite-basalt and intermediate-felsic, calc-alkaline igneous rocks that forms a linear belt along the western margin of the Gindalbie Domain. The greenstones are overlain by late basin sequences, which are dated at younger than ~2673 Ma, are mainly comprised of turbidite with subordinary carbonaceous shale, sandstone, conglomerate, chert and magnetic shale [7,72,84,85,87].

4. Gold Mineralization

4.1. Gold Endowment

In the EGST, gold was first discovered at Coolgardie in 1892 followed by Kalgoorlie in 1893. Since then, more than 20 deposits with metal contents of greater than one million ounces (Moz) of gold (Au) have been discovered in the region. Production to date exceeds 145 Moz Au with the Golden Mile at Kalgoorlie accounting for 65 Moz Au (or ~45%) of the current total [78,88]. Combining past production and unmined resources, which are estimated at greater than 125 Moz Au [1], gives a total endowment of the EGST of well over 270 Moz Au. The study area captures >40.0 Moz Au of the endowment, largely contained in 11 gold deposits and deposit clusters [1] (Figure 1 and Figure 3, Table 2).

4.2. Gold Deposit Styles

The EGST is home to a diverse range of gold deposit styles that are controlled by a variety of different structures and hosted by a variety of rock types of different metamorphic grades from sub-greenschist to amphibolite facies. Despite this diversity, their many common features strongly indicate that the gold deposits in the EGST, and broader Yilgarn Craton, form part of regionally extensive gold mineralizing systems that operated at various times during the late Archean evolution of the host craton [103,104].

According to Tripp et al. [105], gold mineralization occurs as structurally and/or lithologically controlled lodes that typically take the form of quartz ± carbonate veins, breccias and disseminations enveloped by gold-related wall rock alteration. Gold occurs in all rock types, although basalt and dolerite are the most frequent and large granitoid bodies are the least frequent hosts [16,105]. The principal control on gold mineralization was structure, across scales, having controlled both hydrothermal fluid flow as well as the architecture and location of gold deposition [105].

4.3. Gold Depositional Events

The EGST experienced multiple, protracted gold mineralization events during the late Archean, related to both extensional and compressional deformation episodes within a 60 m.y. timeframe between ~2670 and 2630 Ma [5,7,8,9,12,15,105]. Whilst gold was deposited during most EGST deformation events [7], formation of the larger deposits can be attributed to three events (Figure 2) [7,8,15]:

• D₃ (~2665 to 2655 Ma): Formation of EGST crustal architecture during a period of NE-SW-directed extension and associated metamorphic core complex formation. Most of the large crustal-penetrating fault systems had already been established at this time, connecting the upper crust to a metasomatized mantle as indicated by the first arrival of mafic and syenitic granites. The mantle link would have resulted in the addition of significant heat into the base of the crust and aided in mantle-to-crust metal transfer. A prominent D₃ deposit example in the study area is Gwalia (>8.2 Moz Au).
• D₄ (~2655 to 2650 Ma): The extensional crustal architecture established during D₃ was inverted during D₄, a brittle-ductile deformation event that can be divided into an initial phase of ENE-WSW-directed shortening (D_4a) followed by a switch in contraction to WNW-ESE-directed shortening (D_4b). Deformation was accompanied by a switch from high- to low-Ca granite magmatism due to crustal melting. Most importantly, D₄ coincided with the most significant gold mineralization event in the EGST. In particular, strike-slip deformation and reactivation of pre-existing structural heterogeneities during D_4b served as a highly effective fluid focusing mechanism. The most prominent D₄ deposit example is the Golden Mile (>65 Moz Au), located ~35 km S of the study area.
• D₅ (~2650 to 2635 Ma): A further switch in the stress field triggered NE-SW-directed shortening associated with dextral strike-slip, thrusting and low-Ca granite melt emplacement. Gold mineralization was controlled by brittle structures. A prominent D₅ deposit example is Sunrise Dam (>10.3 Moz Au), located ~20 km E of the study area.

5. Data Integration and Interpretation

5.1. Insights from the Enhancement Filtering of Geophysical and Remote Sensing Data

Enhancement filtering of gravity data delivered two critical insights regarding proximity, association and abundance relationships between gravity features and gold occurrences:

• Belt-parallel gravity ridges serve as proxies for greenstone basin centers or zones of thicker mafic-ultramafic material along greenstone keels and root zones (Figure 4). As demonstrated by a confidential gravity data enhancement filter for the entire Australian continent, a large proportion of Yilgarn gold (and nickel) deposits have proximity, association and abundance relationships with gravity ridges. A similar relationship can be observed within the study area where gold occurrences, other than those wholly contained within intrusive rocks, have proximity, association and abundance relationships with gravity ridges observable at multiple scales (Figure 5).
• Gravity edges, on the other hand, represent discontinuities such as fault or shear zones, or significant lithological boundaries. In the case of the study area, the Thunderbox, King of the Hills, Gwalia and Ulysses gold deposits all sit along the same curvilinear belt-parallel gravity edge, which can be traced along strike for >230 km and appears to continue north beyond the study area boundary. Interestingly, this gravity linear only coincides with the Perseverance-Ockerbury-McClure fault system in places whilst there is strong spatial coincidence between it and the basement granite-greenstone contact (Figure 6). Similar, semi-parallel linear features in the study area may represent attractive targets for first pass exploration.

Magnetic edges typically represent fault or shear zones or lithological boundaries. Lithostructural features within the study area that are clearly evident in magnetic edge enhancement filters data are (i) Proterozoic dolerite dykes, (ii) certain intrusive contacts, and (iii) segments of major and subsidiary shear and/or fault zones. In addition, all of the larger gold deposits are located along and coincident with magnetic edges (Figure 7).

Where outcropping at surface or exposed in open cut mines, orogenic gold systems in the study have proximity, association and abundance relationships with (i) domains of high (95th percentile) K/Th values as extracted from radiometric data, and (ii) goethite-clay-iron feature depth index (IFD) ratio anomalies as extracted from remote sensing imagery.

Overall, the application of advanced geophysical and remote sensing data modeling served to map potential ore-controlling structures and ore signatures, thereby greatly enhancing the toolbox for targeting gold in central-eastern Yilgarn Craton. The semi-automated filtering approaches applied to the geophysical and remote sensing data helped mitigate some of the bias and inconsistencies that can occur when data are exclusively interpreted by humans.

5.2. Solid Geology (Lithostructural) Interpretation

The main objective for undertaking a lithostructural interpretation was the generation of a coherent base map for MPM and target generation. Precedence was given to the development of a new structural framework that integrates the results of the geophysical enhancement filtering. Whilst our interpretation incorporated and honored aspects of previous solid geology maps produced by the Geological Survey of Western Australia, pmd*CRC and St Barbara Limited (Table 1), in particular geology, the overall structural framework was rebuilt from the ground up, informed by the new geophysical evidence.

In the study area, form lines (Figure 8a) trace stratigraphy and schistosity in greenstone, compositional layering in granitic and schistosity in gneissic units. Form line orientations are highly variable. However, NNW-SSE orientations dominate, reflecting the principal stratigraphic and tectonic grain of the greenstone sequences. Local variations are mostly due to folding at the kilometer to multi-kilometer scale.

The major, first- and second-order faults within the study area (Figure 8b-d), likely the products of D₃ and subordinate to earlier-formed, N-S-striking, D₁ master faults (cf. Figure 3 in [7]), include the (i) Menzies-Boorara fault system, which separates the Ora Banda and Boorara domains of the Kalgoorlie Terrane, (ii) Ockerburry-McClure-Perseverance fault system, which separates the Kalgoorlie and Kurnalpi terranes, (iii) Melita-Emu fault system, which separates the Gindalbie and Menangina domains of the Kurnalpi Terrane, (iv) Keith-Kilkenny fault system, which separates the Menangina and Murrin domains of the Kurnalpi Terrane, and (v) Celia fault system, which separates the Murrin and Edjudina/Laverton domains of the Kurnalpi Terrane.

The subsidiary third- and fourth-order fault-fracture dataset shows several structural corridors with N-S, E-W, NE-SW and NW-SE orientations dominant (Figure 8e-f). These corridors are oblique, or orthogonal, to the generally NNW-SSE-trending greenstone belts, and may have played an important role in focusing fluid flow and controlling the location of gold mineralization. The latter thesis is supported by the results of a Fry analysis, which identified NNW-SSE, N-S to NNE-SSW, NE-SW and NW-SE directions of maximum gold occurrence alignment (Figure 9).

Lithologically, the study area can be divided into two broad domains, both of which are Archean in age (Figure 10): (i) Felsic to intermediate granite-gneiss, and (ii) greenstone terrain. Intrusive rocks are the most abundant rock type, covering c. 62% (in plan view) of the study area. They include large regional as well as some deformed, gneissic (‘external’ or ‘basement’) granitoid plutons, commonly of monzogranitic composition and often clearly mapped by well-defined gravity lows. Smaller ‘internal’ plutons that intruded the greenstone terrain range in composition from granitic to gabbroic. The greenstones, often spatially coincident with gravity highs, were assigned here to six broad lithological groups:

• Late basin sequences: Polymictic to oligomictic conglomerate and sandstone.
• Sedimentary siliciclastic sequences: Felsic volcaniclastic sandstone, siltstone ± conglomerate, chemical sedimentary and interleaved felsic volcanic rock.
• Felsic to intermediate volcanic sequences: Rhyolite to rhyodacite, andesite, felsic volcaniclastic and siliciclastic rock ± felsic and mafic volcanic rock, interleaved with coeval basalt and dolerite.
• Mafic to ultramafic volcanic sequences: Basalt, komatiite, peridotite, serpentinite, mafic to ultramafic schist (chlorite schist, tremolite schist, fuchsite-andalusite schist) ± coeval dolerite and gabbro.

Figure 8. (continued). Lithostructural (solid geology) interpretation. First- and second-order fault systems draped over (c) gravity (0 to 6400 m differential upward continuation of the Bouguer anomaly grid) and (d) magnetics (1VD of the RTP).

Figure 8. (continued). Lithostructural (solid geology) interpretation. First- and second-order fault systems draped over (c) gravity (0 to 6400 m differential upward continuation of the Bouguer anomaly grid) and (d) magnetics (1VD of the RTP).

Figure 8. (continued). Lithostructural (solid geology) interpretation. (e) Fault and fracture systems. (f) Fault and fracture systems color-coded according to dominant orientation, thereby highlighting several corridors with dominant N-S, E-W, NE-SW and NW-SE structural orientations.

Figure 8. (continued). Lithostructural (solid geology) interpretation. (e) Fault and fracture systems. (f) Fault and fracture systems color-coded according to dominant orientation, thereby highlighting several corridors with dominant N-S, E-W, NE-SW and NW-SE structural orientations.

Figure 9. Fry analysis. (a) Structural framework and gold occurrences used for Fry analysis. (b) Fry plot showing translations of all gold occurrences (i.e., relative center-to-center distances and directions). (c) Rose diagram of all translations illustrating frequency distributions and directions of maximum gold occurrence alignment: NNW-SSE, NNE-SSW, NE-SW and NW-SE.

Figure 9. Fry analysis. (a) Structural framework and gold occurrences used for Fry analysis. (b) Fry plot showing translations of all gold occurrences (i.e., relative center-to-center distances and directions). (c) Rose diagram of all translations illustrating frequency distributions and directions of maximum gold occurrence alignment: NNW-SSE, NNE-SSW, NE-SW and NW-SE.

Within all the above sequences there are narrow, curvilinear units characterized by strong positive magnetic anomalies, typically comprised of banded iron formation (BIF), basalt, dolerite, gabbro, komatiite, peridotite or ultramafic schist.

Last but not least, the study area is traversed by numerous E-W- to ENE-WSW-striking dolerite intrusions of the Paleoproterozoic (~2408 to 2401 Ma) Widgiemooltha dyke swarm [107,108] (Figure 10). The dykes are easily identifiable in magnetic datasets due to their high magnetic contrast and linear nature. Two groups of dykes are evident: (i) 2408 ± 3 Ma dykes that produce positive magnetic anomalies and typically strike c. 075°, and (ii) 2401 ± 1 Ma dykes that produce negative magnetic anomalies and typically strike c. 085°. Their contrasting magnetic polarities are linked to their emplacement during different paleomagnetic regimes. The dykes, which comprise of solidified mantle melts, are interpreted here to have exploited pre-existing crustal heterogeneities and, as such, may mark deep-seated, mantle-tapping structures of Archean age that were re-activated and re-used during Proterozoic times. Interestingly, many of the gold occurrences in the study area are located along or proximal to (i.e., within 5 km of) the dyke swarm corridors, in particular the larger gold deposits.

Figure 10. Lithostructural (solid geology) interpretation—final map product. (a) Simplified version showing undivided granite and greenstone units. (b) Full version showing lithological subdivisions. (c) Legend to Figure 10b.

Figure 10. Lithostructural (solid geology) interpretation—final map product. (a) Simplified version showing undivided granite and greenstone units. (b) Full version showing lithological subdivisions. (c) Legend to Figure 10b.

6. Targeting Model

Generating an efficient targeting model requires careful consideration of several interrelated factors [21,30,109,110]:

• The geologic processes needed to form the targeted mineral deposit type [20,27,48,111,112].
• The mappable expressions of these processes in available geoscience and exploration data [48,113].
• The scales at which these processes have occurred and are mappable [113,114].
• Geological uncertainty and human bias [27,115].

Here we followed the mineral systems approach to targeting [21] with our exploration targeting model structured around the critical processes of orogenic gold deposit formation and their mappable expressions as summarised in Table 3.

7. Mineral Potential Modeling (MPM)

7.1. Statistical Evaluation of Predictor Maps

The predictive performance of each predictor map was statistically assessed using a reference dataset consisting of 1773 gold occurrences, referred to here as mineral deposit locations (MDLs), and 1773 unmineralized sites, referred to here as non-deposit locations (NDLs). The MDLs are comprised of gold deposits, occurrences, and gold mineralized drill intercepts. The NDLs were selected based on the following criteria: They must (i) be spatially randomized; (ii) be situated at a sufficient distance from known prospect sites; and (iii) located outside of the gold permissive tract.

Two statistical approaches were utilized to evaluate the predictor maps: The (i) normalized density index (N_d) [126,127] illustrated in prediction-area plots (Figure A1) and area under the receiver operating characteristic curve (AUC) (Figure A2) [127,128]. Predictor maps that satisfied the statistical thresholds of N_d > 1.00 [124] and AUC > 0.50 [128] were deemed suitable for incorporation into MPM (Table 4).

7.2. Continuously-Weighted Mineral Potential Models

Gold prospectivity in the Leonora district was modeled using three continuous integration techniques: (i) Fuzzy gamma [64], (ii) geometric average [65] and (iii) improved index overlay [66]. Predictor maps were weighted through logistic [63] and small fuzzification functions [129,130], yielding fuzzy membership scores within the interval [0,1] for the proxy variables (Figures A5–A36).

The MPM outcomes for the Leonora district, derived from the application of these integration functions, are depicted in Figures 11a–c.

7.2.1. Continuous Fuzzy Gamma Technique

Continuously-weighted fuzzy spatial predictors can be integrated through the application of fuzzy operators [24]. In the present study, the fuzzy gamma operator that synthesizes the SUM and PRODUCT operators was utilized to achieve fine-tuned calibration of various inputs [129]:

$${\mu }_C=[1-{\displaystyle \prod _{i=1}^n(1-{\mu }_i){]}^{\gamma }\times [{\displaystyle \prod _{i=1}^n{\mu }_i{]}^{1-\gamma }}}$$

(1)

For each cell, ${\mu }_i$ represents the ith input spatial predictor’s fuzzy score, while ${\mu }_C$ denotes the resultant potential score obtained through the combination procedure, with the parameter constrained to the interval $(0\le \gamma \le 1)$. Here, a gamma value of 0.9 was adopted.

7.2.2. Geometric Average Technique

The geometric average constitutes an effective multi-criteria decision-making (MCDM) framework, operationalized for the coherent synthesis of continuously-weighted predictor maps in MPM. Pursuant to the mathematical formulation in Equation 2, the geometric average (G_A) associated with each spatial cell is formally computed as the nth root of the overall product of the respective proxy values [65]:

$$G_A(F_1,F_2,\mathrm{...},F_n)=\sqrt[n]{{\displaystyle \prod _{i=1}^n{F}_i}}=\sqrt[n]{F_1{F}_2\mathrm{...}{F}_n}$$

(2)

For each spatial cell, n denotes the total number of predictor maps, while F_i represents the fuzzy weight corresponding to the ith predictor map. Accordingly, Equation 2 was reformulated to accommodate the specific conditions of the study area as follows:

$$G_{Gold}(F_1,F_2,\mathrm{...},F_{28})=\sqrt[28]{{\displaystyle \prod _{i=1}^{28}{F}_i}}=\sqrt[28]{F_1{F}_2\mathrm{...}{F}_{28}}$$

(3)

In this formulation,$F_1,F_2,\mathrm{...},F_{28}$ denote the fuzzy scores assigned to the evidential values within their respective predictor maps. Subsequent to the computation of $G_{Gold}$ for each unit cell throughout the study area, the resulting values were spatially mapped to develop the geometric average gold potential model.

7.2.3. Improved Index Overlay Technique

The Improved Index Overlay represents an enhanced version of the data-driven index overlay approach applied in both brownfield and greenfield MPM [66]. In contrast to conventional and data-driven index overlay approaches, this technique assigns weights to each predictor map using a modified Shannon’s entropy algorithm. Notably, it eliminates reliance on expert opinion or pre-existing mineral occurrences. The detailed workflow of the modified Shannon’s entropy weighting procedure is outlined in [66]. Following this procedure, the weights of individual predictor maps were objectively determined (Table 5), demonstrating that the density of lithological contacts represents the most influential predictor map, whereas proximity to principal faults exhibits the lowest relative significance.

Upon deriving the weights of the predictor maps through the modified Shannon’s entropy algorithm, the maps were subsequently integrated using the Improved Index Overlay framework as follows:

$$I_{IO}={\displaystyle \frac{{\displaystyle \sum _i^n{W}_{vi}{W}_i}}{{\displaystyle \sum _i^n{W}_i}}}$$

(4)

where for all cells within the study area, $I_{IO}$ denotes the Improved Index Overlay score, while $W_i$ signifies the weight of the ith predictor map, assigned utilizing the modified Shannon entropy algorithm. The variable $W_{vi}$ reflects the continuously allocated cell value of the ith predictor map through a logistic or small fuzzification function. In this study, Equation 4 is formulated as follows:

$$I_{IO}={\displaystyle \frac{W{V}_1{W}_1+W{V}_2{W}_2+\mathrm{...}+W{V}_{28}{W}_{28}}{W_1+W_2+\mathrm{...}+W_{28}}}$$

(5)

Here, $W_1,W_2,\mathrm{...},W_{28}$denote the weights of the predictor maps, determined based on the modified Shannon entropy algorithm. The parameters $W{V}_1,W{V}_2,\mathrm{...},W_{28}$ represent the continuously assigned weights reflecting the evidential value of each individual cell within the respective predictor maps.

7.3. Knowledge-Driven BWM-SAW Technique

Within the spectrum of knowledge-driven techniques to MPM, multi-criteria decision-making (MCDM) approaches are particularly valued for their efficiency [62,132,133,134,135]. In this context, the assignment of weights to predictor maps is generally carried out through either matrix- or comparison-oriented procedures. Among these, the simple additive weighting (SAW: [68]) and best-worst method (BWM: [70]) represent prominent examples of matrix- and comparison-oriented techniques, both of which have been successfully implemented in MPM [62,132,133,134,135].

In this study, a hybrid MCDM framework, namely BWM–SAW [134], was employed in conjunction with the overall performance (O_p) index [136,137,138]. The O_p index is calculated on the basis of an improved P–A plot that incorporates three fundamental criteria: (i) occupied area curve, (ii) NDL prediction rate curve and (iii) MDL prediction rate curve [136]. To achieve this, the BWM technique was implemented to objectively quantify the weights of the decision criteria encapsulated by the predictor maps (Table 4). The O_p index was utilized to determine the weights of the predictor maps, as well as to identify the best and worst predictor maps (Figure A33, Table 6). Subsequently, the others-to-worst (OW) and best-to-others (BO) vectors were established (Tables A1 and A2) in accordance with the O_p values reported in Table 6. Based on these parameters, problem 3 of [134] can be formulated as presented in the following Equation 6:

$$\begin{array}{l}\min \xi \\ \mathrm{s}.\mathrm{t}.\\ \left|W_1-2W_2\right|\le \xi ,\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{all}\hspace{0.17em}j\\ \left|W_1-2W_3\right|\le \xi ,\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{all}\hspace{0.17em}j\\ \left|W_1-2W_4\right|\le \xi ,\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{all}\hspace{0.17em}j\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\\ \left|W_1-9W_{28}\right|\le \xi ,\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{all}\hspace{0.17em}j\\ \left|W_2-5W_{28}\right|\le \xi ,\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{all}\hspace{0.17em}j\\ \left|W_3-5W_{28}\right|\le \xi ,\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{all}\hspace{0.17em}j\\ \left|W_4-5W_{28}\right|\le \xi ,\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{all}\hspace{0.17em}j\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}⋮\\ \left|W_{27}-2W_{28}\right|\le \xi ,\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{all}\hspace{0.17em}j\\ {\displaystyle \sum _j{W}_j=1}\\ W_j\ge 0,\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\mathrm{all}\hspace{0.17em}j\end{array}$$

(6)

Next, optimal weights $(W_B^{\ast },W_1^{\ast },\mathrm{...},W_W^{\ast })$ and ${\xi }^{\ast }$ were determined (Table A3) by solving Equation 6. For $a_{BW}=a_{13}=9,$the obtained consistency index was 5.23 [134] whilst the consistency ratio was $0.009/5.25=0.002,$demonstrating appropriate consistency. After assigning weights to the predictor maps, the SAW approach was employed to rank the alternatives. To achieve this, an initial decision matrix $B_{1201497\hspace{0.17em}\times \hspace{0.17em}28}$was developed that includes 1,201,497 decision alternatives, each linked to an individual cell with a particular coordinate in the corresponding predictor maps, and 28 decision criteria. The alternatives were then ranked via the SAW step-by-step process, previously described by [134]. The BWM-SAW gold potential model for Leonora district is shown in Figure 11d.

7.4. Data-Driven Random Forest (RF) Algorithm

A growing number of recent investigations have underscored the efficacy of RF-based MPM, consistently evidencing its outperformance relative to a broad spectrum of supervised machine learning techniques previously adopted for MPM [31,50,62,139,140,141,142,143]. The superior performance of RF is largely attributable to its inherent ability to alleviate overfitting tendencies and to enhance predictive efficiency through the incorporation of the bootstrap aggregation (bagging) paradigm [67]. RF is an ensemble-oriented machine learning algorithm that capitalizes on a resampling mechanism to generate multiple random subsets of training samples, each subsequently employed in the construction of an unpruned decision tree [67]. Within this framework, a bootstrapping procedure is operationalized, wherein sub-samples with replacement are drawn from the spatial proxy values at mineral deposit and non-deposit locations, collectively designated as labeled data. Decision trees are trained on the in-bag portion of these samples, representing approximately two-thirds of the labeled dataset, while the residual fraction, referred to as out-of-bag (OOB) samples, is exploited to estimate decision tree impurity. As such, RF represents an aggregation of unpruned decision trees, each of which is developed in accordance with a distinct set of so-called existent patterns.

The implementation of RF requires the specification of two principal hyperparameters: (i) the number of decision trees (n) to be generated within the ensemble and (ii) the number of predictor variables (m) considered at each node split [67]. In this study, an n value of 1,000 and m value of 9 were considered based on the procedure outlined by [139]. The relative importance of predictor variables can be quantitatively evaluated through the mean decrease in the Gini impurity index (M_G) and accuracy (M_a). The M_G statistic reflects the degree to which each variable contributes to enhancing node homogeneity and, by extension, the overall discriminative power of the RF model. Conversely, M_a is derived from the computation of out-of-bag (OOB) error, thereby serving as an accuracy-based measure of predictor relevance. In line with [67], higher values of both M_a and M_G indicate greater predictive importance of a given variable. Figure A34 illustrates the significance of the predictor variables in the RF modelling pursuant to M_G and M_a. According to this figure, proximity to felsic to intermediate volcanic rocks is the least significant predictor map in this study while the proximity to domains of juvenile crust, based on the M_a, and the proximity to basement granitoids, based on the M_G, are the most significant ones.

Figure 11e illustrates the RF-based gold potential model generated with the 27 competent fuzzified predictor maps (Figures A5b–A31b), deliberately excluding the predictor map of “proximity to known gold occurrences” to mitigate exploration bias, and a training error curve as depicted in Figure A35. As evidenced by the latter, the modeling error rate exhibits a progressive decline with the incremental growth of decision trees, decreasing from an initial mean squared error of ~0.065 for the first tree to about 0.028 upon the 1,000th iteration.

Figure 11. (continued). Gold potential maps. (e) Random forest.

Figure 11. (continued). Gold potential maps. (e) Random forest.

8. Discussion

8.1. Geology and Structure

One of the key geological observations in this study is that of gold occurrences having proximity, association and abundance relationships with structures that are orthogonal, or at high angles, to the dominant NNW-SSW grain of the terrane- and domain-bounding fault systems and regional schistosity preserved in the greenstone belts.

Firstly, we observed a spatial relationship between gold occurrences, in particular the larger ones, and E-W-striking magnetic edges, representing certain geological features with Proterozoic dolerite dyke corridors the most prominent of these (Figure 7). A similar relationship was noted by [144] for the Yilgarn Craton as a whole in which most of the major gold districts exhibit “intriguing spatial associations” with Proterozoic dykes. This association was interpreted by [144] in terms of the Proterozoic dykes occupying much older planes of weakness that formed prior to or at the time of Archean gold mineralization. We agree with this interpretation and believe that Proterozoic dyke swarms, which can be easily traced in magnetic ± gravity data, map important yet cryptic basement structures, some or segments thereof played a critical role in localizing gold mineralization. This belief was confirmed by spatial statistical analyses undertaken as part of this study, returning AUC values of 0.54 and 0.52 as well as N_d values of 1.08 and 1.04 for predictor maps representing Proterozoic dolerite dyke proximity and density, respectively. According to published constraints [126,128], any predictors with AUC > 0.50 and N_d > 1.00 are statistically sound, in this case confirming a valid spatial relationship between gold occurrences and Proterozoic dolerite dykes, or more accurately the structures that host them.

Secondly, we observed N-S, E-W, NE-SW and NW-SE structural trends, represented by subsidiary, late tectonic (likely of D_4/5 timing), low-displacement faults and fractures, that are oblique, or orthogonal, to the dominant NNW-SSW grain. Fry analysis indicated similar directions of maximum gold occurrence alignment, implying that structures of this nature may have played an important control on the location of gold mineralization. Some of these structures are coincident with NE-SW- and NW-SE-striking low-displacement faults and fractures interpreted by [145], who referred to them as Perkins Discontinuities. The latter term describes planar boundaries that mark an abrupt termination of mineralization continuity, typically observed at the hand specimen to deposit scales. According to [145], the nature of the Perkins Discontinuities identified at the Aphrodite deposit and surrounding area suggests that gold mineralization occurred after the formation of these discontinuities, resulting in asymmetric gold mineralization across these low-displacement faults and fractures, which themselves are typically unmineralized. A lithostructural interpretation by [146] of the Bardoc area, done in greater detail than the solid geology interpretation completed as part of this study, showed that Aphrodite lies within a ~5 km-wide, NE-SW-striking corridor defined by a series of fractures and low-displacement (typically <100 to 300 m) faults with apparent dextral kinematics (Figure 12). These structures, as well as less common NW-SE-striking ones, cut obliquely across the dominant NNW-SSE trend of the Bardoc Tectonic Zone. Deposit-scale studies completed at Aphrodite reported cross structures of broadly similar orientation, including NE-SW-striking and steep to moderately NW-dipping faults, the NW-SE-striking and moderately NE-dipping Epsilon orebody [93] and abundant yet poorly mineralized, ENE-WSW (070°) and WNW-ESE (130°)-striking quartz-carbonate ± pyrite veins [93,147]. We believe that the intersection of the dominant NNW-SSE and subordinate NE-SW (± 22.5°) and NW-SE (± 22.5°) structural trends created a permeability network that played a major role in localizing gold mineralization, not only at Aphrodite but also in other places where similar structural interaction and complexity exist. More work would be required, and at a variety of scales and locations, to better understand how cross structures influenced gold mineralization and the geometry thereof.

8.2. Geophysics and Remote Sensing

Magnetic and gravity data are the most commonly used and available geophysical datasets for regional targeting of orogenic gold mineralization in the Yilgarn Craton. Greenstones often show significant contrast relative to the surrounding granite in both magnetic and gravity data. Detailed magnetic data typically map the fabric of greenstone belts. Fabric-parallel features apparent in magnetic data tend to be fault or shear zones or stratigraphic boundaries. By and large, ‘cross-fabric’ features are widespread with the most common ones typically representing cross-cutting fault-fracture systems and mafic dikes. Gravity data over the Yilgarn Craton are generally too coarse in resolution for mapping detailed features. However, they are useful for mapping greenstone distribution and relative thicknesses. They also serve to highlight major cross structures. The main goal of the semi-automated geophysical filtering undertaken as part of this study was to highlight geological features in a relatively unbiased and consistent way. One way this was done is by extracting features in particular orientations (Figure 5, Figure 6 and Figure 7). This can be useful because different orientations often are related to different styles or periods of deformation. The phase congruency algorithm extracts features in eight different orientations (i.e., 0°, 22.5°, 45°, 67.5°, 90°, etc.). These can then be combined to highlight structures in certain orientation ranges. For example, the fabric parallel features in Figures 5a–7a are generated by combining the 0° and 167.5° orientation so that features with a NNW-SSE to N-S trend are highlighted. In the case of the Yilgarn Craton, these features tend to be belt-parallel fault/shear zones or stratigraphic contacts. Gold mineralized structures are typically parallel to the fabric of the host greenstone belt. Cross structures often impact grade continuity and distribution. Extracting these different types of structures separately makes targeting more effective. The phase congruency method is also relatively amplitude independent, which makes it possible to highlight subtle features in the data. This is important because orogenic gold systems are not always associated with strongly magnetic bodies. Rather, the geophysical signatures of orogenic gold systems are highly variable due to a wide range of host rocks, structural settings, mineralization styles, ore minerals and alteration assemblages.

Direct detection of unweathered orogenic gold systems with Sentinel-2 data is difficult. Pyrite is one of the most common alteration phases, but it cannot typically be detected with satellite remote sensing data. However, the Sentinel-2 satellite can detect the weathering products of pyrite such as goethite and jarosite [45,148]. Weathering of pyrite generates acidic conditions that lead to enhanced generation of clays during the weathering of silicate minerals [149]. The Yilgarn Craton has experienced significant weathering and oxidation [150], so any outcropping gold mineralization should be weathered to an assemblage containing ferric-Fe-bearing phases and clays. Areas with coincident highs in goethite, clay and IFD indexes were found to correlate well with known mineralization in the Leonora area. Any single remote sensing product will have many false positive anomalies. We helped limit this issue by using a CMY ternary combination of the goethite, clay and IFD indexes. Our approach allowed us to highlight all areas with coincident highs in these indexes accurately and efficiently. Extracting the anomalies using an automated approach helps mitigate issues with interpreter bias or inconsistencies. The main limitation of this approach with the Sentinel-2 data is that it can only detect outcropping mineralization. Even very thin cover would block any signal from underlying rocks.

8.3. Mineral Potential Modeling (MPM)

8.3.1. Predictor Map Performance

The significance of predictor maps in MPM can be conceptualized from two perspectives. First, during the pre-modeling stage, attention is directed toward the identification and selection of the most competent predictor maps (Table 4), which are subsequently incorporated into the integration framework for developing the potential model. Second, in the post-modeling stage, predictor maps are systematically ranked according to the degree to which they have contributed to the prediction of potential values (Figure A34). From the pre-modeling perspective, the careful selection of well-performing predictor maps plays a critical role in minimizing systemic uncertainty and ensuring that the modeling process is grounded on robust and meaningful spatial variables. From the post-modeling perspective, predictor maps that exhibit limited or negligible predictive capacity may be excluded from further consideration, thereby enabling the modeling process to be refined and repeated with a reduced but more effective set of inputs. Such an iterative approach not only facilitates the optimization of the modeling framework but also allows for the comparative evaluation of successive mineral potential models, providing insights into the sensitivity of model performance to the choice of predictor maps.

It is strongly recommended to use a combination of statistical analysis tools to identify effective predictor maps and to quantify their contribution in predicting potential values. This approach enables the comparison of predictor maps based on different performance metrics, thereby allowing for more accurate decision-making regarding their differentiation, elimination or retention. In the former perspective, for instance, whilst both the “proximity to fold hinges” (predictor M1) and “proximity to flanks of granitoid bodies” (predictor M2) predictor maps achieved N_d scores of 1.440, their AUC scores set them apart in that predictor M1 scored an AUC value that is ~1.24 times higher than that of predictor M1. Thus, considering only the N_d index would not have allowed us to differentiate between the two predictor maps; however, incorporating the AUC index highlighted the superiority of predictor M1 over predictor M2. As another instance, as illustrated in Figure A34, the predictor map of “proximity to areas of demagnetization” is ranked 21st according to M_a, but 9th according to M_G. To achieve a more precise comparison, it is advisable to incorporate additional performance metrics, such as entropy or information gain.

8.3.2. Multi-Technique Approach

The robustness of MPM is highly contingent upon the quality of the underlying genetic model and translation of the key genetic processes into mappable targeting criteria, since the conceptual model sits at the core of the MPM process and stochastic uncertainty mostly arises from inappropriate selection of predictor maps. Whilst these issues can be minimized by way of mineral systems modeling, systemic and stochastic uncertainties can be further mitigated by adopting a multi-technique approach to MPM, which (i) ensures optimal use of the available empirical and conceptual information, (ii) compares, contrasts and cross-validates MPM results and (iii) better constrains exploration targets by integrating the results from diverse numerical models [31,114]. As evidenced in [31,50,62], a multi-technique approach to MPM not only facilitates the generation of more robust exploration targets but also provides insights unattainable through any single modeling technique, while also contributing to the development and calibration of novel tools and techniques.

Beyond the aforementioned advantages, adopting a multi-technique approach coupled with the development of diverse gold potential models (Figure 11) also supports appraisal of model performance. To achieve this, we applied the improved P-A plot procedure of [136]. Pursuant to Figure A36 and the performance statistics reported in Table 7, the RF algorithm to MPM delivered the best-performing gold potential model for the study area (overall performance, O_p = 0.52) (Figure 11e). It is imperative to note that the potential model developed by the RF algorithm demonstrates different patterns compared to the ones developed by the other techniques (Figure 11). This difference, improving the reliability of prediction in this study, is an intrinsic property of the RF algorithm, which constructs an ensemble of regression trees through bootstrap sampling and derives the final prediction by averaging their outputs, thereby reducing variance and enhancing robustness.

8.3.3. Predictor Map Sensitivity

Several of the predictor maps used in this study have potential to introduce bias into the MPM given they (i) present ‘self-fulfilling prophecies’ (e.g., proximity to gold occurrences) by enhancing prospectivity in areas around the known gold deposits, (ii) have possible interdependencies with other predictors (e.g., proximity to greenstone belts or proximity to the flanks of granitoids) given spatial overlap, (iii) are ineffective in areas of cover (e.g., proximity to remotely sensed alteration systems or proximity to domains of high K/Th values) because certain signals do not penetrate the regolith, and/or (iv) do not uniformly and objectively represent the study area as a whole (e.g., proximity to fold hinges) because of inconsistent mapping scales data resolutions.

Given that such biases can have serious consequences [114], we set out to investigate model sensitivity to certain predictor maps. More specifically, we employed the best-performing RF approach to generate two additional gold potential maps for the study area, altering the number and types of input predictor maps: Alternative RF model 1, which excludes predictor maps (a) proximity to gold occurrences, (b) greenstone belts, (c) flanks of granitoids, (d) fold hinges, (e) remotely sensed alteration systems and (f) domains of high K/Th values, and alternative RF model 2, which excludes predictors (a), (b) and (c). Both alternatives recorded excellent overall performance (O_p) values in the range from 0.49 (alternative RF model 1) to 0.54 (alternative RF model 2), respectively. These values are close to or, in the case of alternative RF model 2, even exceed the O_p value of the original RF model of 0.52. The results are not only similar statistically but also visually with little discernible difference between the original and alternative RF-generated gold prospectivity maps (Figure 13).

8.3.4. Using Fractal Thresholding for Target Generation

The concentration-area (C-A: [151]) fractal technique, a proven and widely used method for classifying continuous spatial data [28,66,127,152,153,154], was applied here to prioritize areas of gold potential derived from the superior RF model. Two thresholds (Figure A37), marking the most significant changes in slope and reflecting higher potential scores as well as data populations more strongly associated with mineralization, were determined at potential scores of 0.830 (2nd order target areas) and 0.958 (1st order target areas). These two thresholds were then used to generate a ternary-class exploration targeting map, in which the modeling results were categorized into domains of high, moderate, and low gold potential (Figure 14). The delineated target areas are spatially correlated with the known gold occurrence locations and key predictors of the underlying genetic and target models for orogenic gold mineral systems. As illustrated in Figure 14, the prioritized high potential, 1st order target area represents an order of magnitude (i.e., >13 times) reduction of the search space, the hallmark of a well-performing, practically useful targeting technique [114].

In this study, apart from identifying high-priority targets, the results obtained from MPM served three additional critical purposes: (i) providing context and constraints for area prioritization of future gold exploration activities, (ii) identifying the most prospective ground held by competitors, thereby flagging areas of interest with regard to potential corporate investment, acquisition or joint venture deals and (iii) providing context and constraints for the monitoring of mineral tenement status and acquisition of new, vacant ground.

8.4. Targeting

8.4.1. Role of MPM as a Targeting Tool

Since its inception in the late 1960s and early 1970s [113], MPM has come a long way both from a method development perspective and as an effective targeting tool. According to [49], whilst methodological and technological improvements have been addressed comprehensively and published widely, demonstrations of practical applications of MPM have been very limited in the literature apart from the showcasing of mineral potential maps. Examples of MPM case studies that have gone a step further in that the modeling results served as input for additional analysis and/or particular attention was paid to the development of comprehensive targeting models and quality targets include [28,49,59,155,156,157,158,159].

As summarized by [31], the discovery in 2020 of the Gonneville Ni-Cu-PGE deposit in Western Australia represents the first world-class discovery that can, at least partly, be attributed to MPM, validating the effectiveness of MPM as a targeting tool and presenting proof-of-concept for this technology. Gonneville was a greenfields discovery by Chalice Mining Limited [160], made on the back of a Ni-Cu-PGE prospectivity map produced by Geoscience Australia [161] that also unlocked an entire new search space for Ni-Cu-PGE deposits along the >1000 km-long, underexplored western margin of the Archean Yilgarn Craton.

Regardless of these abilities, MPM generated prospectivity maps do not guarantee discovery success and, thus, should not be treated as ‘treasure maps’. Rather, they should serve as starting points for ‘treasure hunts’. As such, MPM is best treated as just another tool in the ‘exploration toolbox’, more specifically, a decision-support tool for delineating, ranking, and prioritizing exploration targets based on the underlying modeled prospectivity. Also important to understand is that prospectivity maps present a snapshot in time, reflecting the conceptual understanding of the targeted mineral deposit type as well as the availability and quality of the supporting data at the time of modeling. New information and/or new or better-quality data will invariably change the results, typically for the better [28].

8.4.2. Assessment of MPM Target Areas

The RF model, the best-performing of the MPM approaches employed in this study, correctly identified the locations of all the large (>1 Moz) gold deposits and deposit clusters in the study area like Paddington, Gwalia, Mt Morgans, Tarmoola/King of the Hills, Thunderbox, Apollo Hill, Aphrodite, Ulysses and Menzies (Figure 13 and Figure 14). It also correctly identified most of the smaller gold occurrences, such that the ratio of deposits located within anomalies to the area covered by anomalies [50,140,162] is 83.93% to 7.69% in the high gold potential domain (red polygons in Figure 14), demonstrating the excellent performance of the RF potential model.

Given these results, the modeled areas of gold potential that have no known gold occurrences but similar mappable ingredients as areas that are gold mineralized clearly merit further investigation, particularly where underexplored. One way to assess whether a prospective area may be underexplored or adequately tested is by examining the available, previously collected exploration data. The most useful data for an initial assessment are drilling data given that drilling provides a definitive test of the subsurface. This is particularly true for reverse circulation and diamond core holes. Areas of great(er) drilling density typically denote more intensive and detailed exploration activity that would have led the previous explorers to develop a reasonable understanding of the subsurface geology of the tested area. In places where exploration was intense and drilling dense and relatively deep (>100 to 200 m, or more) it is unlikely that a gold deposit would have been missed if present.

Open-file drilling data for the study area were sourced largely from the Geological Survey of Western Australia (GSWA), augmented by additional confidential drilling data provided by St Barbara Limited (SBM) (Table 1). The GSWA database comprises of >231,000 records whilst the SBM dataset covers >77,000 records (Table 8). Some overlap exists between these two databases. In addition, the GSWA records are open-file only, which means that any as-of-yet confidential drill holes are missing. Despite these shortcomings, the GSWA and SBM drill hole databases provide a reasonably complete picture of previous drilling completed within the study area. Whilst, at first glance, the greenstone belts in the study area appear to be well tested (Figure 15b), this picture changes when drill collars other than those relating to reverse circulation and diamond core drilling are removed. When only plotting the collars of the deeper reverse circulation and diamond core drill holes it becomes evident that large tracts of greenstone are untested at depth (Figure 15c).

Our review of the previous drilling (Figure 15) identified six broad areas of high gold potential that have relatively low drill hole densities, in particular when only considering reverse circulation and diamond core drill holes with hole lengths of 100 m, or more. These areas feature many of the targeting parameters outlined in Table 4, and, therefore, present attractive targets.

A manual analysis [59,163] of these areas identified eight more discrete targets within the broader prioritized MPM target areas (Table 9, Figure 16). Apart from having received less drilling compared to neighboring areas, all but one of the manual targets are located along NNW-SSE-striking, ‘800 m-scale’ gravity ridges. This ridge filter maps the thicker portions of greenstone belt that often correlate well with significant gold deposits. As illustrated in Figure 5, many of the large gold deposits occur on or adjacent to 800 m-scale ridges, in particular those that trend NNW-SSE.

The top-two ranked manual targets are particularly exciting:

• Target #1: Covers a cluster of poorly tested intrusions of the McAuliffe Well Syenite that intruded mafic volcanic-dominant greenstone sequences, abut the first-order order Keith-Kilkenny fault system and are cut by E-W- to NW-SW-striking dolerite dykes. Shallow (maximum hole depth of 45 m) rotary air blast drilling by Saracen Gold Mines Proprietary Limited returned intercepts up to 1.00 m @ 12.28 g/t Au from 9.00 m (hole YER143), 1.00 m @ 6.94 g/t Au from 10 m (hole YER158) and 2.00 m @ 2.20 g/t Au (hole YER159), defining surficial, saprolite-hosted gold mineralization over an area of 400 by 300 m at the Dingo prospect. The adjacent Bull Terrier prospect returned up to 16.00 m # 2.59 g/t Au from 105.00 m (hole YBD-2) and 16.00 m @ 1.52 g/t Au from 60.00 m, including 1.00 m @ 211.70 g/t Au from 65.00 m (hole YRC-63) [164,165]. No deeper or broader, more systematic exploration drilling appears to have been completed across the McAuliffe Well Syenite intrusive cluster. Moreover, the cluster is held by various exploration companies. The disjointed ownership, which appears to have prevailed over decades, likely played a role in preventing a more holistic exploration approach.

• Target #2: Covers prospective mafic and felsic volcanic-dominant greenstone sequences, comprising chemically reactive banded iron formations (BIF) and syenite intrusions, along strike from the Mt Morgans gold production center. Previous drilling returned drilling intersections of up to 6.70 m @ 13.15 g/t Au from 95.00 m (hole MRC036), 5.90 m @ 7.24 g/t from 79.00 m (hole MRC003) and 2.90 m @ 5.41 m from 112 m (hole MRC028) from and defined a small (>150 koz Au @ 1.40 g/t Au) resource at the Korong-Waihi prospect. The latter remains open along strike and at depth with no prior drilling below a vertical depth of 150 m [166,167]. Similar to target #1, ownership of target #2 is disjointed, currently preventing a broader, more holistic approach to exploration.

9. Summary and Conclusions

This study of Archean orogenic gold systems in the central-eastern Yilgarn Craton illustrates that a geophysically-informed lithostructural framework, when combined with a mineral systems-based targeting model and a suite of complementary prospectivity algorithms, can provide a powerful tool for gold exploration.

Enhancement filtering of gravity, magnetic, radiometric and remote sensing data revealed cryptic structures and alteration systems that correlate with known gold deposits but extend into underexplored ground.

The solid geology map generated in this study underpinned a multi-technique mineral potential modelling (MPM) workflow comprising five complementary prospectivity techniques spanning knowledge-driven, data-driven and continuous methods. The best performing random forest (RF) correctly predicted the majority of known gold occurrences while generating a greatly reduced, prioritized search space for new gold targets.

The key strengths of our approach are its reproducibility, explicit incorporation of often cryptic deep-seated architecture and cross-structure networks, and ability to integrate and test multiple weighting schemes to minimize bias. The results confirm that integrating a detailed geophysical interpretation with mineral systems-based prospectivity modelling is an effective way to identify concealed gold corridors and reduce exploration risk. The methodology is readily transferable to other commodities and geological provinces worldwide.