Predictive targeting in Australian orogenic-gold systems at the deposit to district scale using numerical modelling

Abstract

3D numerical models of coupled deformation and fluid flow provide a useful tool for exploration in orogenic-gold systems. Numerical modelling of ore-forming processes can lead to a reduction in targeting and detection risk, thus improving the value proposition of mineral exploration. Hydrothermal mineralisation arises from a complex interplay of deformation, fluid flow, conductive and advective heat transport, solute transport and chemical reactions. Coupled simulation of all of these processes represents a significant computational challenge that cannot be solved within the time-scale of a mineral exploration program. However, the problem can be simplified by identifying a subset of processes representing the first-order controls on mineralisation at the scale of interest. For most orogenic-gold systems, it is argued that the first-order controls on mineralisation at the camp to deposit scale are deformation-induced dilation, fluid flow and fluid focusing. Hence, numerical models of coupled deformation and fluid flow can provide a quantitative insight into the localisation of ore-forming fluids in this type of system. In two case studies, known deposits were modelled in order to determine the critical deformation and fluid-flow-related factors controlling the localisation of mineralisation in these systems. The quantitative results from the forward models were then used as a basis for constructing predictive models that were applied to regional targeting, prospect ranking and selecting the choice of detection methods. Both case studies show that numerical modelling is capable of reproducing the distribution of known anomalism, and that it can predict anomalies that were not expected or accounted for by purely empirical analysis.

Key Words:

• Australia
• deformation
• fluid-flow
• Kundana
• mineralisation
• numerical modelling
• orogenic gold
• Stawell

Introduction

This paper aims to demonstrate the value of integrating 3D numerical modelling of coupled deformation and fluid flow into the task of orogenic-gold-system exploration. Two case studies illustrate the positive impact that numerical modelling of geological processes can have on the mineral exploration value proposition, through its potential to improve prospect, target and detection technique selection. We discuss the validity of, and the approach to, applying 3D deformation and fluid-flow modelling to orogenic-gold systems, briefly review previous applications of numerical modelling to ore-systems analysis, and provide two case studies illustrating the application of numerical modelling to two different exploration problems.

Mineral exploration value proposition

The task of mineral exploration is fundamentally about maximising the probability of discovery while minimising the discovery cost. One way of separating out the different components of risk inherent in this statement is to break it down into the four major questions which need to be addressed by a company assessing investment in a region (Roberts & Eadie 2000; Penney et al. 2004): (i) what is the probability that the target orebody is present in the target area; (ii) if it is present, what is the probability that it can be found using the planned exploration strategy; (iii) if it can be found, what is the probability that the ground will be accessible to the company; and (iv) what is the probability that an action by government will render the orebody uneconomic? The focus in this paper is on the first two questions, which can be summarised as targeting and detection risk, respectively, together with a consideration of the cost involved in using detection methodologies.

Targeting risk can be defined as the probability that an area selected for application of some detection technology (e.g. drilling, geophysics, geochemical sampling, etc.) will not contain the target ore or, at larger scales, the target orebody. Detection risk can be defined as the probability that the ore cannot be found in the volume of rock of interest by the detection methodology chosen. Detection risk can generally be reduced by using more and more expensive methods (e.g. intensive drilling).

There is a view that targeting risk is reduced at smaller scales and therefore that the focus of ‘predictive discovery’ should be at terrane scale and above. While the premise is clearly true, there is an interplay between targeting and detection risk which means that targeting decisions remain highly important even at quite small scales.

At larger scales, the discovery challenge becomes greater. Stepping out to the district scale, the other focus of this paper, a major reason for discovery failure in many smaller companies is ‘gamblers ruin’: that is, what happens when explorers are obliged to cease exploration because they have exhausted their funds before discovering an orebody that is present on their tenement holding (Singer & Kouda 1999)? At this scale, both targeting and detection risks are in play: targeting risk—because the explorer may be testing the targets in the wrong order, and therefore gamblers ruin might result, and detection risk—because the explorer is using a detection technique that does not detect the orebody reliably, or which generates many false positives that are expensive to follow up, also leading to gambler's ruin.

Thus, if numerical modelling is to add value to mineral discovery rates, it must either improve prospect selection or help to identify the most cost-effective detection technique. The case studies presented in this paper demonstrate how this can be achieved.

Rationale for deformation—fluid-flow modelling of orogenic-gold systems

The general problem of hydrothermal ore formation involves coupled interactions between deformation, fluid flow, conductive and advective heat transport, solute transport, and chemical reactions (Figure 1). Simulation of this fully coupled problem presents a significant computational challenge. Currently, there is no code that is capable of simulating the fully coupled problem over geologically relevant time-scales in a system with realistic physical and geochemical complexity, within a manageable computation time. A loose coupling between mechanical behaviour and chemically reactive transport in multiphase, non-isothermal systems (Rutqvist et al. 2002) has been applied to numerical simulations of carbon dioxide sequestration in aquifers and to disposal of nuclear waste in unsaturated fractured porous media. This coupled approach has also been used in a study of the mechanical and hydrothermal effects of the injection of deep-sourced fluids into a hot, shallow hydrothermal system (Todesco et al. 2004), but has not been applied to any mineral exploration-related research.

Figure 1 Feedback between processes involved in the formation of hydrothermal ore deposits.

To date, there is no published work describing the successful application of fully-coupled 3D mechanical/thermal/fluid-flow/reactive transport (geochemical) modelling to mineralisation processes. However, identifying the first-order controls on mineralisation for a given system, and limiting the simulation to this subset of the controlling processes, can provide useful insights so long as the chosen subset is appropriate to the system under consideration. The most appropriate subset will depend on the mineralisation style, the tectonic setting of the system, and the scale at which it is modelled.

Garven and co-workers (Garven & Freeze 1984a, b; Ge & Garven 1992; Raffensperger & Garven 1995; Appold & Garven 1999, 2000; Garven et al. 2001; Simms & Garven 2004) chose the coupled thermal – fluid-flow – geochemical subset when modelling topographically driven and buoyancy-driven systems such as the Mississippi Valley-type base-metal deposits (Leach et al. 2001; Huston et al. 2006), and the Athabasca Basin unconformity-related uranium deposits of northern Saskatchewan, Canada (Hoeve & Sibbald 1978; Ruzicka 1995). Xu et al. (2001) modelled supergene copper deposits, Schardt et al. (2005) the VHMS deposits of the Pilbara Craton in Western Australia, Yang et al. (2004) the HYC stratiform Pb – Zn deposit of the Northern Territory and Matthai et al. (2004) the Mt Isa copper deposits in Queensland. Kühn et al. (2006) and Kühn & Gessner (2006) also simulated the Mt Isa copper deposits where they applied coupled 3D thermal, fluid-flow and reactive transport modelling. All of these studies used coupled modelling of fluid flow and heat transport and/or geochemistry, as these were considered to be the first-order controls on the mineralisation in these systems. With the exception of Kühn et al. (2006) and Kühn & Gessner (2006), geometries were simplified to 2D representations of the geology. Deformation was not represented in any of these models.

This paper focuses on the application of 3D numerical modelling to the simulation of processes responsible for the formation of gold-only deposits in orogenic-gold systems such as those described by Groves et al. (1998) and Goldfarb et al. (2001), and reviewed by Goldfarb et al. (2005) and Robert et al. (2005). For these systems, the initial premise is that the structural- or deformation-related controls on dilation, fluid-flow and fluid focusing are the critical first-order controls on the localisation of mineralising fluids, and therefore the localisation of gold mineralisation, at the deposit to camp scale. While not all favourable structural sites will host mineralisation, it is argued that without the influence of deformation-enhanced fluid focusing, a significant orogenic-gold deposit cannot form. It is acknowledged that magmatic or thermally driven fluid-flow can play a significant role in determining the location of such deposits at the crustal or craton scale, and that geochemical controls are critical at the scale of the lode or vein. However, the weight of evidence indicates that the primary requirement for a significant orogenic-gold deposit is a structurally controlled focus for mineralising fluid flow (Groves et al. 2000; Robert et al. 2005; Bierlein et al. 2006).

The case studies presented in this paper involve deposit- to camp-scale numerical models of orogenic-gold systems in which the gold is hosted predominantly in veins or other structures formed by shear and/or tensile failure. The approach to modelling these systems is based on the following assumptions about orogenic-gold systems at the deposit to camp scale: (i) sites of maximum fluid flux and/or fluid focusing are assumed to be good proxies for vein-hosted gold mineralisation, where gold dominantly co-precipitates with quartz; (ii) deformation-induced dilation and permeability enhancement are assumed to be the dominant controls on the fluid-flow regime; (iii) the contribution of thermal buoyancy to the driving force for fluid flow is assumed to be negligible relative to the effect of deformation on fluid flow; and (iv) temperature variations are insignificant and do not provide insight into the location of mineralisation on this scale. Given these assumptions, the decision was made to model coupled deformation and fluid flow, and not include heat transport, thermal buoyancy, solute transport or chemical reactions in the simulations. It is noted that fluid focusing towards dilatant sites may be a good proxy for mineralisation because rapid dilation may lead to rapid fluctuations in fluid pressure (Sibson et al. 1988), potentially resulting in fluid unmixing (Groves & Foster 1991). Dilatant sites may also provide locations for fluid mixing, which was proposed by Hagemann et al. (1994) and Walshe et al. (2003) to be a potential catalyst for gold precipitation. However, we do not set out to show that any of these mechanisms is the dominant control on mineralisation.

Numerical modelling of coupled deformation and fluid flow has been applied to several types of hydrothermal ore systems, including the tin deposits of the Malage orefield in China (Jiang et al. 1997); the Mary Kathleen U – REE deposit (Oliver et al. 1999, 2001) and the Cu and Pb – Zn – Ag deposits of Queensland (Oliver et al. 2001; Ord et al. 2002; Zhang et al. 2006b); the Yilgarn gold province (Sorjonen-Ward et al. 2002) and Hamersley iron-ore system (McLellan et al. 2004) of Western Australia; gold mineralisation at Bendigo, Victoria (Schaubs & Zhao 2002), at Stawell, Victoria (Schaubs et al. 2006), and in the Hodgkinson Province of Queensland (Vos et al. 2007); Cu – Au deposits of Papua New Guinea (Gow et al. 2002); the Outokumpu Cu – Zn – Co deposit, Finland (Ord & Sorjonen-Ward 2003; Zhang et al. 2006a); and the deposits of the Shuikoushan district, China (Zhang et al. 2007). The basis for much of this work lies in understanding and simulating the localisation of deformation (Hobbs & Ord 1989; Hobbs et al. 1990) and its effect on fluid flow in the upper crust, where deformation is assumed to be represented by a pressure-dependent elastic – plastic constitutive behaviour (Ord 1991).

One alternative modelling approach that has been applied to orogenic-gold systems is known as stress-transfer modelling. This is a mechanics-only modelling method which does not simulate fluid flow. It involves calculating the instantaneous (static) stress change arising from a specified slip event on a pre-defined fault, assuming poro-elastic behaviour, and then using this stress change to identify areas that have been brought closer to failure. It has been shown that these areas correspond closely to the location of known orogenic-gold deposits associated with some crustal-scale strike-slip fault systems (Cox & Ruming, 2004; Micklethwaite & Cox 2004, 2006). This relationship may be explained by an increase in permeability in areas brought closer to failure (Sheldon & Micklethwaite 2007). The static nature of stress-transfer modelling makes it relatively quick and easy to use. Hence, it can be a useful exploration tool in cases where mineralisation can be related to specific, small-magnitude slip events on pre-defined surfaces, and it has been shown to be predictive in such systems (Cox & Ruming, 2004; Micklethwaite & Cox 2004, 2006). However, the usefulness of this technique is limited by the static nature of the stress-transfer calculation, and by the assumption of poro-elastic behaviour. Stress-transfer modelling does not represent time-dependent effects associated with fluid flow and moving boundaries. Furthermore, the stress calculation becomes invalid if the stress state moves beyond the yield envelope, which is likely to occur for anything other than a single, small-magnitude slip event. Finally, stress-transfer modelling cannot be used to explore the effect of a specified far-field stress on a system with complex architecture and varying mechanical and/or fluid-flow properties. These issues can be resolved by modelling elastic – plastic constitutive behaviour coupled with fluid flow, as in the case studies presented in this paper.

Stress mapping is another modelling approach that has been applied with some success to orogenic-gold systems (Ojala et al. 1993 Holyland & Ojala 1997, Mair et al. 2000). This technique uses elastic or elasto-plastic mechanical simulations to predict locations of low minimum stress (σ₃) and low mean stress, and thus the loci of dilation. Stress mapping has yielded some useful insights into the distribution of stress and its impact on dilation and mineralisation in structurally controlled ore deposits, including the Sunrise Dam (Mair et al. 2000) and Granny Smith (Holyland & Ojala 1997) gold deposits in the Yilgarn. However, the usefulness of this technique is limited by the fact that fluid flow is not explicitly modelled; it is only inferred to occur in response to dilation. Thus, there is no quantification of flow rates or volumes, or of changes in effective stress associated with fluid flow.

Theoretical basis and numerical method

The theory, governing equations and numerical method underlying the simulations of coupled deformation and fluid flow described in this paper are the same as those employed in a number of previous studies (Ord & Oliver 1997; Oliver et al. 1999; Gow et al. 2002; Ord et al. 2002; Schaubs & Zhao 2002; Sorjonen-Ward et al. 2002; McLellan et al. 2004; Schaubs et al. 2006; Zhang et al. 2006a, b, 2007). Here, an overview is provided without going into details of the equations or numerical method.

Deformation is assumed to be governed by an elastic – plastic constitutive model, comprising linear elasticity combined with the Mohr – Coulomb yield function and non-associated plasticity. The theory underlying this constitutive model has been summarised by Vermeer & de Borst (1984) and Mandl (1988). The constitutive model is defined in terms of the elastic moduli, yield functions for shear and tensile failure, and corresponding potential functions defining the direction of plastic (irreversible) strain associated with yielding. The yield functions are defined by the cohesion, friction angle and tensile strength. A material conforming to this constitutive model behaves elastically, typically reducing in volume under compressional deformation, until the stress state reaches the yield envelope, at which point it begins to undergo shear or tensile failure and develops plastic (irreversible) strain. Shear failure is associated with dilatancy governed by the dilation angle, representing the effect of initially close-packed grains moving apart as they slide over one another (Vermeer & de Borst 1984; Ord & Oliver 1997). This dilation typically drives fluid flow from elastic and/or compacting regions of the model towards the regions of dilation. Shear strain can become localised to form features similar to geological shear zones, even in an initially homogeneous material. This localisation is a consequence of non-associated plasticity in the Mohr – Coulomb constitutive model (Ord 1991; Vermeer & de Borst 1984). Deformation experiments suggest that the cohesion, friction angle and dilation angle vary with strain and confining pressure (Ord 1991), but it would be impractical to quantify such variations for every modelled lithology at the pressure – temperature conditions relevant to the system being modelled: hence, constant values are used for each rock type in these simulations. This is a reasonable approximation given the uncertainty in rock properties on the metre to kilometre scale and at the PT conditions relevant to the system being modelled.

Fluid flow is assumed to be governed by Darcy's law with isotropic permeability, such that the direction and magnitude of the fluid flux are governed by the hydraulic head (i.e. the fluid-pressure gradient minus the hydrostatic component) and by spatial variations in permeability (Phillips 1991). Fluid properties such as compressibility and viscosity are assigned constant values representative of pure water at Earth surface conditions. As with the mechanical properties, this is a useful first approximation in the absence of detailed information concerning the composition and properties of pore fluid in the system of interest. Deformation and fluid flow are coupled in the following ways: (i) fluid pressure increases with compaction and decreases with dilation; (ii) fluid pressure influences the effective stress, which in turn influences deformation (e.g. an increase in fluid pressure may result in shear or tensile failure); (iii) deformation leads to changes in hydraulic head as the fluid moves with the deforming rock; and (iv) in some models, permeability is increased during shear and tensile failure.

The case studies described in this paper use the Lagrangian finite difference code, FLAC3D to model coupled deformation and fluid flow in 3D (Cundall & Board 1988; Itasca 2002). Simulations are defined in terms of: (i) their architecture or geometry; (ii) the mechanical and fluid properties assigned to the mesh elements (e.g. cohesion, friction angle, permeability); and (iii) the boundary conditions considered most relevant to the geological scenario represented by the simulation (e.g. boundary velocities and pore pressures representing a given tectonic regime). Each of these attributes can be varied to test their effect on the system. Values for the mechanical properties may be guided by geotechnical data, if available, or estimated based on the relative strengths inferred for the pre-mineralisation rock mass and published data for similar rock types. The models are initialised with supra-hydrostatic pore pressures, consistent with interpretations of deformation features such as tension veins in many orogenic-gold systems (Groves et al. 1998; Sibson et al. 1988; Mair et al. 2000; Goldfarb et al. 2005). Subsequent deformation modifies the distribution of stress and fluid pressure, bringing some regions of the model closer to the failure envelope and moving others further from it.

Faults or shear zones may be defined in the initial model architecture as weak zones of finite width, and additional shear zones may develop during deformation due to spontaneous localisation of shear strain. Computational limitations inhibit accurate representation of structures that are very small or narrow relative to the overall size of the model: hence, modelled fault zones tend to be wider than the real structures that they represent. This approximate representation of narrow structures may impact on the absolute values of stress, strain and fluid flux predicted by a given model, but is unlikely to have a significant impact on the overall pattern of behaviour. Hence, interpretation should focus on the relative differences, or qualitative aspects, of the results and comparisons between model scenarios, rather than on absolute values of any parameter from a specific model.

The pattern of deformation and fluid flow that develops within a model reflects the combined effects of the boundary conditions and evolving conditions in the model. For example, uniform upward flow representing a fluid source at depth can be imposed by specifying a uniform fluid flux or elevated fluid pressure at the lower boundary of a model. This uniform flow field may be perturbed by localised dilation associated with shear or tensile failure within the model. The impact of localised dilation on fluid flow is a key factor in the case studies presented in this paper.

Episodic slip associated with fluctuations in fluid pressure, permeability and mechanical properties is a common feature of fault behaviour in the brittle upper crust (Cox et al. 2001). The models presented in this paper do not generally predict such behaviour: instead, faults or shear zones in the models tend to fail continuously once the stress state has reached the yield envelope, for as long as the stress state remains on the yield envelope. This continuous failure mode may be thought of as a time-averaged approximation of multiple, discrete slip events, or aseismic creep. However, episodic behaviour can arise in these models as a consequence of the interplay between permeability, fluid pressure and effective stress.

The scalar product of the fluid flux with the pressure gradient, u•▿p, is employed as an analysis tool in one of the case studies presented in this paper. This parameter can be used to highlight regions of enhanced fluid flux in the direction of decreasing fluid pressure, such as regions of fluid flow towards dilatant sites. The value of u•▿p becomes increasingly negative with increasing fluid flux in the direction of decreasing fluid pressure. Negative values of u•▿p may indicate precipitation of pressure-sensitive minerals such as quartz, which is known to co-precipitate with gold in orogenic-gold systems, assuming pressure is the dominant control on solubility at the scale of interest (Phillips 1991). In practice, it is difficult to quantify deposition rates of specific minerals without carrying out a full reactive transport simulation, taking account of variations in temperature and bulk chemical composition as well as fluid pressure. The use of u•▿p as an analysis tool does not imply or require that fluid pressure is the dominant control on gold mineralisation; rather, it is used in conjunction with other indicators to aid in visualisation of the fluid flow regime (note that u•▿p is a scalar parameter that can be contoured), and in particular to identify locations where fluid may be focused and flowing rapidly towards dilatant sites.

Modelling approach

The application of coupled deformation and fluid-flow numerical modelling to mineral exploration problems typically follows a two-stage approach. The two stages are known as the model validation phase and the predictive modelling phase. Initial simulations are designed to test the conceptual geological/mineralisation model(s) for a known deposit in a region of interest: this is the model validation phase. These simulations aim to reproduce, or forward model, the observed distribution of mineralisation, strain or alteration within the known system in order to constrain the critical parameters controlling this distribution. In orogenic gold and other structurally controlled ore systems, such parameters may include, but are not limited to: far-field stress orientation or deformation regime; rheological or permeability contrasts between lithological units within the system; geometric attributes of the system such as the orientation of faults or stratigraphic units; and the degree of fluid overpressure. These parameters are regularly discussed and described in the existing literature on orogenic gold and structurally controlled ore systems. In some cases, constraining the importance of a particular parameter, for example the importance of rheological contrasts associated with bodies that have anomalous geophysical properties, can lead to the preferential acquisition of data to directly detect prospective locations.

The model validation phase not only yields constraints for these key parameters, but also gives an indication of the sensitivity of the system to variations in these parameters. Another potential outcome of the model validation phase is the identification of previously unknown or unexpected anomalies, or the identification of preferential structural host orientations, within the ore system. This can yield either a better understanding of the grade distribution within the deposit, or identify new structural orientations or locations requiring a revised detection method (such as a different drill direction) in order to optimise detection and definition.

The quantitative results of the model validation phase constrain the input parameters for generic or specific predictive models, which aim to predict the location of anomalies outside the area of the known deposit. The two case studies below illustrate both the model validation and predictive stages of modelling.

Examples

The two examples illustrate the application of numerical modelling to orogenic-gold systems at Stawell in Victoria, and at Kundana in Western Australia. The intention is not to present exhaustive case studies of these gold systems, but to illustrate briefly the application of two different numerical modelling approaches to two different deposit types, and to show how both the model validation phase and the predictive phase impacted on the exploration process.

Case study 1: Stawell, Victoria

The information from which this summary was drawn is described in substantially more detail by Schaubs et al. (2006). This study was based on earlier work by researchers from the University of Melbourne (Miller et al. 2001; Miller & Wilson 2004) and geologists from Leviathan Resources Limited.

Geology

The Stawell Corridor and the Magdala mine are located in a reworked Cambrian orogenic zone located between the Lachlan and Delamerian Orogens (Miller et al. 2005) (Figure 2a), and known as the Stawell Zone. In this area, basalt domes with associated gold mineralisation occur in a region ∼15 km wide, bounded on the west by the east-dipping Moyston Fault and on the east by the west-dipping Coongee Fault (Figure 2b).

Figure 2 (a) Simplified geological map of southeastern Australia focusing on the Western Lachlan Fold Belt (modified after VandenBerg et al. 2000). (b) Geological map of area surrounding Stawell (modified after Schaubs et al. 2006).

The Stawell Corridor extends to the north-northwest under Murray Basin sediments (Figure 2b). These increase gradually in depth northward, and for this reason, the area was essentially unexplored prior to the work of MPI Mines and Leviathan Resources. Nonetheless, the basalt domes persist under cover along with substantial potential for gold orebodies similar to those at Stawell itself.

The domes are made up of massive to pillowed tholeiitic basalt (Squire & Wilson 2005) and are overlain by the Stawell Facies (Dugdale et al. 2006). The latter comprises calcareous and quartzose sandstones, chert and mudstone, some of which is sulfidic and is interpreted to have been altered prior to mineralisation (Watchorn & Wilson 1989). East of the Coongee Fault, the St. Arnaud beds (Figure 2b) are made up of less deformed and less metamorphosed sandstone and shale.

At the Stawell mine, gold mineralisation in the Magdala orebody occurs primarily on the west flank of the Magdala Dome, a northwest-trending doubly-plunging basalt body ∼4000 m long (Figure 3a, b). On the flanks and top of the dome, there are a number of basalt lobes that resemble parasitic folds (Figure 4), but are thought to represent primary flow lobes formed during the extrusion of the basalt (Squire & Wilson 2005). The Stawell Facies overlies the basalt dome (Miller & Wilson 2002). Metasedimentary rocks surrounding the Magdala Dome are divided into two units: the Albion and Leviathan Formations (Squire & Wilson 2005). The former is more pelitic and the latter more psammitic (Miller & Wilson 2002).

Figure 3 Location of mineralisation at the Dukes Nose and Golden Gift relative to the simplified reconstructed Magdala Basalt Dome (pre-D₆ movement on the South Fault). (a) View towards the southeast of the entire dome. (b) View from above the dome. (c) View towards the southeast of mineralisation in the immediate area of the Dukes Nose Basalt lobe (modified after Schaubs et al. 2006).

Figure 4 Cross-section through the Magdala antiform at section 320. Note location of Dukes Nose and Extended Basalt lobes on the west flank of the dome (modified after Schaubs et al. 2006).

The deformation history at the Magdala mine is well documented (Miller & Wilson 2002, 2004). Three early ductile deformation events pre-dated the gold mineralisation and produced a variably developed layer-parallel schistosity, upright folds with a strong axial-planar fabric and a differentiated crenulation cleavage and refolding, respectively.

Gold mineralisation was coeval with two subsequent brittle deformation events. The first (D_4a) resulted in the development of northeast-striking reverse faults due to dominantly east-northeast – west-southwest-directed contraction which is interpreted to have rotated to east – west (D_4b). The second event (D₅) marked a switch to a sinistral shearing environment characterised by tension gashes near the basalt that formed during northwest – southeast-oriented shortening. The northwest-striking and northeast-dipping South Fault (Figure 4) is interpreted to have been active during D₆ deformation and offset the Magdala Antiform (Miller et al. 2001; Miller & Wilson 2004).

Gold occurs in a number of structural settings and locations in the Magdala ore system, including the Central and Golden Gift Lodes. This case history focuses on mineralisation which occurs within the Stawell Facies and associated structures on the margin of the basalt; little mineralisation is present in the basalt itself (Robinson et al. 2006). In this structural and stratigraphic position, the bulk of past and present resources are located in the Central and Golden Gift Lodes. The Central Lode occurs on the west flank of the dome at the contact between the Stawell Facies and the Albion Formation at the Dukes Nose (Figure 3c). The Golden Gift Lodes occur on the offset portion of the Magdala Antiform below the South Fault (Figure 3a) (Miller et al. 2001; Miller & Wilson 2004).

Simulation

Model validation

The model validation phase tested the simple proposition that the location of gold mineralisation was most influenced by the shape of the basalt dome and its rheological contrast with the adjacent Stawell Facies. It was recognised that the chemical composition of the Stawell Facies probably exerted an additional control on ore formation (Dugdale et al. 2006). The Stawell Facies is not equally gold-mineralised everywhere, so there must have been another overriding control on ore localisation. The numerical modelling experiments described herein sought to explain ore formation in terms of focused fluid flow and dilation, and to test the proposition that the first-order control on ore localisation was the orientation and magnitude of the pore fluid pressure gradient (p) with respect to the local fluid-flow vector.

The model validation phase required relatively limited experimentation with the various parameters which might influence dilation and fluid flow because: (i) the relative strengths of the rock units were derived from earlier geotechnical testing of rocks at the Magdala Mine (these measurements were consistent with inferred relative strengths of the rock units at the time of mineralisation, although it is important to note that this may not always be the case as subsequent alteration, deformation or weathering can alter absolute rock strength characteristics post-mineralisation); (ii) the stress history of the Stawell area had been extremely well characterised by earlier work (Miller et al. 2001; Miller & Wilson 2004); (iii) the 3D geometry of the Magdala basalt dome in the vicinity of the mine workings was well known through drilling; and (iv) the post-ore fault movement history was sufficiently well understood to enable reconstruction of the shape of the basalt dome at the time of mineralisation. Despite this, some 50 models were run to test the known range of different stress orientations on the various basalt dome shapes (Schaubs et al. 2006).

A table of the rheological properties employed in the models described in the following section is provided as Table 1. The tensile strength values used are consistent with typical values of 1 – 10 MPa for rock (Etheridge 1983; Sibson & Scott 1998). Permeability remained constant during deformation. The model was fully saturated, and no external fluid sources or fluxes were applied. Pore pressure values were initialised to a pore fluid factor of 0.5, halfway between hydrostatic and lithostatic. The presence of quartz tension veins at Magdala suggests that overpressured conditions existed locally and transiently, but that the region was not subjected to supra-lithostatic pressures continually. The depth of burial suggests that purely hydrostatic fluid pressure conditions were unlikely at the timing of mineralisation. All models were deformed to 5% shortening as the structural observations of Miller & Wilson (2002) suggested that the mineralising events at Stawell were associated with small strain, which is consistent with previous conclusions that low syn-mineralisation displacements are common in orogenic-gold deposits (Ojala et al. 1993).

Table 1  Stawell model: rock properties (Schaubs et al. 2006).

	Basalt	Albion Formation	Magdala Facies
Density (kg m^−3)	2950	2760	2790
Bulk modulus (GPa)	67	45	61
Shear modulus (GPa)	45	20	22
Cohesion (MPa)	40	30	25
Tensile strength (MPa)	4	3	3
Friction angle (°)	30	25	20
Dilation angle (°)	2	3	4
Permeability (m^2)	1 × 10^−16	1 × 10^−15	5 × 10^−12

Models were run at two scales, covering just one basalt lobe (the Dukes Nose, the source of much of the historic gold production) and the scale of the entire reconstructed basalt dome (Figure 3a). Results from two of the simulations are shown in Figures 5 (small lobe scale) and 6 (dome scale). The overall distribution of known ore-grade gold mineralisation at both the mine and dome scales correlates moderately well with maximum simulated fluid-flow rates for the D_4a and D_4b stress orientations (Figures 5, 6). Also, simulations at the shoot scale undertaken by Robinson et al. (2006) provided an even closer match to ore-grade distribution than the dome and lode-scale models.

Figure 5 (a, b) Accumulated positive volume strain (dilation) in models of the Dukes Nose after 5% shortening. Only the highest ranges of values (black patches) are shown and are as follows: (a) 0.05 – 0.086; (b) 0.05 – 0.086. Values below 0.05 are transparent. (c, d) Instantaneous fluid-flow velocities in models of the Dukes Nose at 5% shortening. Only highest values (represented by arrows) are shown and are as follows: (c) 4.5 × 10⁻⁷ to 4.97 × 10⁻⁷ m/s; (d) 3.5 × 10⁻⁷ to 4.98 × 10⁻⁷ m/s. Values below 3.5 × 10⁻⁷ m/s are transparent. Note abbreviations (e.g. E – W) for shortening directions.

Figure 6 Instantaneous fluid-flow velocities in models of the Magdala Dome at 5% shortening. Only highest values (represented by arrows) are shown and are as follows: (a) 8.5 × 10⁻⁸ to 1.12 ×10⁻⁷ m/s; (b) 8.5 × 10⁻⁸ to 1.08 × 10⁻⁷ m/s. Left and right columns represent two different views of the same models. Values below 8.5 × 10⁻⁸ m/s are transparent. Note abbreviations (e.g. E – W) for shortening directions.

Maximum simulated fluid-flow rates are particularly favoured where there is a transition from contraction to dilation. This is exhibited in the above models where there is a significant change in the dip of the basalt on the western flank, at the northern end, near lobe crests, and in the vicinity of the basalt lobes, such as those in the Dukes Nose area. Stawell Facies rocks between the lobe and the main dome contacts contract in the simulations, while the overlying rocks dilate, thus creating the required pore-fluid pressure gradient. The fact that this phenomenon will occur is rather intuitive: however, the exact positions of the regions of maximum fluid flux are controlled by subtle variations in basalt shape and are therefore difficult to predict without the use of numerical modelling. Regions of maximum fluid flow in these models do not occur at the very top of the dome, but rather occur further down the flanks; below areas of maximum dilation and low minimum principal stress (σ₃). This differs to the results of stress mapping carried out by Mair et al. (2000) who presented models where low minimum principal stress occurs at the hinges of folds and crests of domes. Both the style and the location of mineralisation with respect to the dome at Stawell differ from the chevron-fold hosted deposits in Victoria (e.g. Bendigo) and those in the Howley district modelled by Mair et al. (2000). In these latter deposits, mineralisation is vein-hosted and occurs at the crests of anticlines which are areas of low minimum principal stress as predicted by both the coupled deformation fluid-flow models presented here and the stress mapping technique.

It is worth noting that these simulations do not reproduce grade distribution perfectly, and the detailed match to ore grade shells becomes poorer as the models become larger. This reflects the need to simplify the geology as the model volume grows. Nonetheless, at each scale, the simulation results provide a 3D map of prospectivity that: (i) encompasses the mineralised zone; (ii) extends beyond the well-characterised (i.e. drilled) area; and (iii) is significantly smaller than the total possible target area (Figures 5, 6). The value of this approach becomes particularly apparent at the dome scale, where a small amount of shallow drilling would have been sufficient to identify the simplified shape of the body, and the simulation could then have indicated the best areas to focus drilling.

Predictive phase: kewell prospect

Once the Stawell model phase had yielded satisfactory results, attention turned to several prospects north of the mine beneath Murray Basin cover, including Kewell and Wildwood (Figure 2). The full results of this work are described in Rawling et al. (2006) and Schaubs et al. (2006). This paper refers to one model only (Figure 7), from the Kewell Prospect.

Figure 7 Fluid-flow velocities in models of the Kewell Dome after 5% shortening. Only highest values (represented by arrows) are shown: (a) 1.0 × 10⁻⁷ to 1.24 × 10⁻⁷ m/s; (b) 8.5 × 10⁻⁸ to 9.84 × 10⁻⁸ m/s. Left and right columns represent two different views of the same models. Values below the low range of velocities are transparent. Note abbreviations (e.g. E – W) for shortening directions. Circles signify where regions of high volume strain, high fluid-flow rate and gold anomalies from air-core drilling coincide.

The geology of the Kewell Prospect is broadly similar to the Magdala Dome area and consists of a basalt dome overlain by Stawell Facies, enclosed by metasedimentary rocks and overlain by post-ore cover sediments ∼100 m thick. Exploration activities prior to the modelling had determined the approximate shape of the Kewell basalt body. Thus, air-core drilling had established the outline of the basalt body immediately below the overlying sediments, and geophysical modelling, based both on gravity and magnetic data, had provided an interpretation of the dips of the basalt margin and hence the plunge of the dome axis at both ends of the body. The shape of the uneroded body was then inferred by projecting upwards from the basalt outline below cover.

The stress history and rheological parameters derived from the model validation phase at Magdala were applied at Kewell. The distribution of simulated fluid-flow vectors suggested that the most prospective areas for gold mineralisation in the uneroded portion of the dome were likely to be at the southern flanks of the dome (Figure 7) and its northern tip. Neither of these areas in the dome margin was highlighted as being prospective at Magdala (Figure 6) in the validation stage, despite the fact that the shape of both domes is broadly similar. In this respect, the results were somewhat counterintuitive.

Exploration completed prior to completion of the predictive simulations also included two diamond drillholes on the western flank of the dome, designed to test air-core geochemical anomalies (Rawling et al. 2006). Neither of these drillholes intersected ore grade and width gold mineralisation, although the second, more southern, hole did encounter more gold anomalism. In retrospect, these results were consistent with the fluid-flow simulation outcomes. However, in the absence of the simulation results, it was not obvious where the next drillholes should be sited.

Subsequent drilling based on the Kewell simulations generated ore grade and width intersections, with the first hole in the program including 4.1 m at 12.6 g/t Au (Dugdale 2004).

Discussion

The Stawell District study demonstrates the potential for obtaining predictive value once a successful model validation step has been achieved. At Kewell, the predictive step yielded results which were not immediately obvious to the mineral explorer and allowed the exploration team to focus an expensive diamond-drilling program on a relatively confined area with a high probability of containing the target gold mineralisation. Thus, in the terms described above, improved prospect quality was obtained through the application of numerical modelling. In addition, the risk of failing to intersect ore-grade gold mineralisation in the next few drillholes was reduced, an important point where the drilling budget is limited, which was the case with the Kewell drilling program (Dugdale 2004).

Case study 2: Kundana, Western Australia

Geology

The Kundana Mineral Field is located ∼20 km northwest of Kalgoorlie, Western Australia. Its gold endowment is substantial, with historic production and resources to 2004 reported to be 217 t of Au (Tripp 2004). The field comprises numerous gold deposits and occurrences which are primarily controlled by localised dilational veining associated with competency contrasts across key stratigraphic contacts (Tripp 2004). The aim of the numerical modelling project was to explore how variations in lithological contact and fault geometries in the Kundana field might influence or control the distribution of mineralisation.

The regional and local geological setting of the Kundana Mineral Field was described by Hadlow (1990) and Tripp (2004), while the ore deposits and mineralisation styles were described by Lea (1998) and Tripp (2004). The following summary is derived primarily from these papers.

The Kundana Mineral Field consists of a series of gold deposits which are generally close to the interpreted position of the Zuleika Shear Zone. The sequence consists of Late Archean intermediate and mafic volcanic rocks, dolerites and sedimentary rocks including graphitic black shale, polymictic conglomerate and felsic sedimentary rocks (Figure 8). The best gold hosts in the field are relatively thin and highly strained graphitic black shale units and andesite – dolerite contacts (Tripp 2004). High-grade mineralisation (up to 90 g/t) is hosted predominantly as coarse gold in thin, planar, laminated quartz veins which dip steeply (55 – 85°) towards the west and are commonly located on stratigraphic contacts. Significant gold deposits in the Kundana field are depicted in Figure 8, and are described in detail by Lea (1998) and Tripp (2004).

Figure 8 Kundana goldfield geology (from Tripp 2004).

The deformation history discussed by Tripp (2004) suggests that a period of coaxial compression or transpression prevailed from the early folding event through to the time at which ore formation is inferred to have taken place. Tripp (2004) proposed progressive overprinting and coaxial reactivation of structures including shear fabrics, lithological contacts/dilation sites and ore-host veins, and suggested that strain partitioning associated with competency contrasts and fabric intersections was instrumental in localising Au-bearing fluids into dilatant veins.

The recognition of a close association of high-grade mineralisation with dilatant stratigraphic contacts and Zuleika-parallel structures and high-strain zones led exploration geologists to target these localities. However, there are several stratigraphic contacts of interest, multiple Zuleika-type structures, and an unresolved association with north- to northeast-trending strike-slip cross faults, within the prospective ∼60 km² of the Kundana district (Figure 9). This resulted in >100 potential targets, with traditional empirical prospectivity analysis methods unable to prioritise these targets with any confidence.

Figure 9 Kundana fault map (courtesy G. Tripp, Barrick Gold).

It was proposed that a generic deformation – fluid-flow modelling exercise be undertaken in an attempt to understand why some fault segments and lithological contacts were mineralised, while others were not. It was hoped that testing a range of fault/contact and cross-fault orientations and a range of far-field stresses would highlight which geometries were most likely to dilate during the ore-forming event. This could then provide a basis for ranking the large number of potential exploration targets within the Kundana field, and potentially predicting previously unrecognised ore trends.

Simulation

The simulations were to be based on the black shale-hosted laminated-vein systems of the Kundana North and South Pits. These systems exhibit a marked competency contrast between a thin, weak, foliated black shale sandwiched between competent basalt in the hangingwall to the west and andesites in the footwall to the east. The key question was one of geometry. The aim was to use numerical modelling to reproduce the observed anomalies associated with mineralised shale and cross-fault geometries at the North and South Pits. The parameters constrained during the forward-modelling phase would then be used to predict sheared shale and cross-fault intersection geometries and orientations which are favourable for dilation throughout the Kundana field.

The variables that were explored in the simulations included: (i) orientation of the stress field at the time of mineralisation; (ii) orientation of the thin sheared shale units (described as F1 below), which were interpreted to have behaved like weak shear zones at the time of mineralisation and were modelled as such; and (iii) orientation of the interpreted strike-slip cross faults (described as F2 below), which were assumed to be present at the time of mineralisation. Figure 10 depicts the model mesh, which was generated using a template (Potma et al. 2004) that facilitates independent variation of F1 and F2 dip and dip-direction as well as far-field stress orientation with respect to the fault geometry.

Figure 10 Template-based geometry and mesh used for the Kundana numerical models involving a parameter search varying F1, F2 and σ₁ orientations. The mesh-building template consists of inner and outer meshes. The model is run on the entire mesh, but the results of interest are only derived from the inner mesh, thereby reducing boundary effects in the model results. The orientation of the inner mesh can be varied with respect to the outer mesh, which enables the modeller to examine the effects of changes in the stress field with respect to a single geological geometry. Orientation, curvature and properties of individual model domains and fault segments are parameterised, enabling rapid mesh generation for a wide range of geometric variants.

As in case study 1, the simulation results were analysed for anomalies in shear and volumetric strain (dilation), failure state and fluid flux. In addition, the scalar product of the fluid flux with the pressure gradient, u•▿p, was used to highlight regions of enhanced fluid flux in the direction of decreasing fluid pressure, such as regions of fluid flow towards dilatant sites.

The range of geometric variables tested in the initial parameter sensitivity study are shown in Table 2 and Figure 11. The key outcomes are summarised below.

Figure 11 Range of F1, F2 and σ₁ orientations applied during the Kundana numerical modelling geometric parameter search.

Table 2  Kundana models: fault geometry parameter search variables.

Primary parameter variables	Far-field stress (σ1)	Sheared shale strike (F1)	Sheared shale dip (F1)	Cross-fault strike (F2)	Cross-fault dip (F2)
Placer Dome—preferred value (best guess)	070°	320°	60°	015°	Vertical
Variable range	060 – 090°	320 – 350°	55 – 85°	000 – 030°	−70 to −110°
Sample interval	10°	10°	10°	10°	10°

Results

Simulations of the weak southwest-dipping foliated black shales (F1) and subvertical north- to northeast-striking cross-faults (F2) successfully reproduced the distribution of strain observed in the shale-hosted deposits. The range of simulated fault intersection geometries (F1, F2) and far-field stress orientations confirmed that variations in these parameters, within the range in the Kundana field, produced variations in shear strain, dilation, fluid-flow rates and anomalous values of u•▿p (95th percentile) which were consistent with the known distribution of strain, alteration, veins and mineralisation. Some representative examples of these 3D model results are presented in Figures 12 and 13. Table 3, ranking the anomalism generated from each of the simulations, highlights that the most critical controls on mineralisation potential are the dip of both the F1 and F2 structures, the strike of the F1 structures and the orientation of σ₁ with respect to the structural architecture, or the overall fault architecture with respect to σ₁.

Figure 12 (a) Impact of varying the principal compressive stress direction on u•▿p anomalism. A 1° change in σ₁ (from 080 to 090°) resulting in a significant u•▿p anomaly. (b) F1 dip impact on dilation within the F1 plane. Note: the scalar values represented by coloured cells have a bottom cut applied at ∼95th percentile (i.e. only the top 5% of values are plotted). The warmer colours represent the extreme anomalism, cool colours still represent significant anomalism, and no scalar colour indicates that anomalism is either less significant or not present. Absolute values are not defined here, as it is the relative difference between models that is the key indicator.

Figure 13 (a) F1 dip impact on u•▿p. Gentle (55°) F1 fault dip resulting in increased u•▿p anomalism within the F2 plane, through a non-intuitive complex feedback mechanism between dilation, fluid-flow velocities and local pore-pressure gradients. (b) Complex feedback resulting from F2 fault dip and dip direction variations creating u•▿p anomalies within the F1 plane. Note: the scalar values represented by coloured cells have a bottom cut applied at ∼95th percentile (i.e. only the top 5% of values are plotted). Colours as for Figure 12.

Table 3  Kundana model results: prospectivity matrix.

Max dilation	u•▿p anomaly	Max dilation	u•▿p anomaly	Max dilation	u•▿p anomaly	Max dilation	u•▿p anomaly	Max dilation	u•▿p anomaly
σ1	σ1	F1 strike	F1 strike	F1 dip	F1 dip	F2 strike	F2 strike	F2 dip direction	F2 dip direction
70	90	140	140	65	65	15	15	90	90
80	80	140	140	65	65	15	15	90	90
90	70	140	140	65	65	15	15	90	90
60	60	140	140	65	65	15	15	90	90
70	70	140	150	65	65	15	15	90	90
70	70	150	140	65	65	15	15	90	90
70	70	160	160	65	65	15	15	90	90
70	70	170	170	65	65	15	15	90	90
70	70	140	140	55	55	15	15	90	90
70	70	140	140	65	65	15	15	90	90
70	70	140	140	75	75	15	15	90	90
70	70	140	140	85	85	15	15	90	90
70	70	140	140	65	65	10	0	90	90
70	70	140	140	65	65	20	10	90	90
70	70	140	140	65	65	30	30	90	90
70	70	140	140	65	65	0	20	90	90
70	70	140	140	65	65	15	15	70 – 285	70 – 105
70	70	140	140	65	65	15	15	80 – 285	80 – 105
70	70	140	140	65	65	15	15	70 – 105	70 – 285
70	70	140	140	65	65	15	15	80 – 105	80 – 285

The orientation of the F1 – F2 architecture with respect to σ₁ has a significant impact on u•▿p values (Figure 12a) and a minor impact on dilation. Of the tested geometries, the tighter the acute intersection angle between σ₁ and F1 strike, the greater the u•▿p and dilation anomalism.

Dilation was noticeably increased on F1 structures where the strike of these structures trended west of north (e.g. an F1 strike of 140° was significantly more dilatant than 170°: Table 3). Similar relationships between mineralisation and fault trend have previously been reported for mineralisation along the Boulder – Lefroy Shear Zone by Weinberg et al. (2004). Models with F1 strikes west of 140° were not tested (at this stage) because F1 structures with this orientation were not known in the Kundana field. However, simulations did predict that if more westerly trending structures were present at the time of mineralisation, they should be highly prospective for Kundana-style vein-hosted gold mineralisation.

The more gently dipping F1 structures were significantly more dilatant (Figure 12b) as one would expect, due to their more preferred orientation for reverse reactivation. However, an unexpected outcome was that F1 dip variation produced significant u•▿p anomaly variation in the F2 plane (Figure 13a). More gently dipping F2 structures also exhibit increased dilation potential in the F2 plane. The dip of F2 structures did not significantly influence u•▿p anomalism. However, F2 dip did influence the nature and distribution of dilation in the F1 plane (Figure 13b), another complex feedback that was not anticipated prior to modelling.

Follow-up predictive simulations indicated that modelling may be able to predict, quantify or simulate other significant mechanisms or behaviours linked to mineralisation within the Kundana field. Some useful outcomes of these simulations are summarised below.

Simulating the presence of a competent unit, such as a gabbro, in the hangingwall of F1 structures resulted in an increased probability of tensile veining and/or brecciation in these units at the time of the mineralising event. It is likely that these veins would be at a high angle to the typical contact-parallel Kundana laminated shear veins because, as tension veins, their orientation will be purely controlled by the orientation of σ₁ and σ₃ with either subhorizontal or approximately east – west- to northeast – southwest-striking vein orientations. This style of mineralisation is present in the immediate footwall sequence of the Rubicon deposit. Reducing the thickness of the F1 shale units in the model, either along strike or down dip, resulted in increased potential for dilation in the thinner portion of the F1 shales (Figure 14a, b), but the u•▿p anomaly was best developed in the thicker portion of the F1 structure.

Figure 14 (a, b) Impact of F1 shale/fault thickness changes on dilation, both along strike (a) and down dip (b). Note that where the F1 structure thins, the dilation anomaly is significantly enhanced. (c) Simulation result analogous to the White Foil deposit, a plunging shoot of weakly distributed dilation to the west of the main Zuleika trend proximal to the south side of a major strike-slip cross-fault. Such an anomaly would be expected to produce a more disseminated or distributed mineralisation style than the conventional narrow high-grade lodes on the Zuleika trend. Colours as for Figure 12.

Some simulations predicted broad regions of distributed dilation and associated u•▿p anomalies in the host rocks surrounding the F1 and F2 structures. One such model resulted in a broad southwest-plunging shoot of distributed dilation near the cross-fault (F2) southwest of its intersection with the F1 structures (Figure 14c). Placer Dome geologists recognised similarities between this model and the relative location, geometry and disseminated mineralisation style of the White Foil deposit. While this single model was not able to predict where other potentially favourable locations for like mineralisation might occur, it was interesting that it unexpectedly reproduced the large-scale morphology and location, with respect to the fault architecture, of the White Foil deposit.

Outcomes

The simulation results indicated which of the known fault geometries within the field were most prospective for Kundana-style high-grade laminated vein formation. These results were validated against known deposits where geometric attributes compared well. Follow-up modelling predicted secondary mineralisation styles and locations within the system, such as discordant tension-vein arrays in competent F1 hangingwall, and potentially footwall, rock packages, and disseminated lower grade mineralisation in the host-rock sequence. These simulations also indicated whether high-grade vein-hosted or disseminated wall-rock alteration/sulfidation reactions are likely to predominate in certain structural positions.

The most significant outcome from the geometric parameter search was the realisation that, if more westerly trending structures were present at the time of mineralisation, they should be highly prospective for Kundana-style vein-hosted gold mineralisation. Placer Dome had previously considered northwesterly to westerly trending faults off the main Zuleika ore trend to be unprospective. However, reinterpretation of magnetic data around an unsuccessful drilling program, which targeted an interpreted Zuleika Shear Zone/stratigraphy-parallel soil anomalies in the south of the Kundana field, resulted in the interpretation of a west-northwest-trending structure. This structure was subsequently drilled, with two diamond drillholes oriented normal to the structure in contrast with the previous drillholes which were oriented at right angles to the north-northwest Zuleika trend. Both holes intersected a significant west-northwest-trending steeply-dipping vein of a similar morphology to the mineralised Kundana veins. Assays of the split core from the two holes returned 0.46 m @ 6.94 g/t Au in one hole and 0.35 m @ 0.35 g/t Au from the second hole. However, the second low-grade interval exhibited visible gold within the remaining split core: thus, the assay is interpreted as an anomalously low-grade attributable to the nuggetty nature of these vein-hosted deposits. Importantly, the first of these two intercepts was higher grade than all of the previous 34 holes drilled within a 500 m-radius of these two holes prior to the numerical modelling exercise. Thus, it was demonstrated that deformation and fluid-flow modelling could be used to identify the optimum target and drill orientation, resulting in ore-grade intercepts, in an area that had previously been extensively drilled without success.

Summary and conclusions

This paper focused on demonstrating the value of applying 3D numerical modelling of coupled deformation and fluid flow to the process of exploring for orogenic-gold deposits.

Hydrothermal ore systems, in general, and orogenic-gold systems in particular, arise from a complex interplay of many factors. The existence of these systems depends on interactions between deformation, fluid flow, conductive and advective heat transport, solute transport and chemical reactions, which in turn are influenced by the tectonic regime and the spatial distribution of lithological units of varying rheology and permeability within the crust (Phillips et al. 1996; Groves et al. 1998, 2000; Weinberg et al. 2004; Goldfarb et al. 2005). Simulation of this complex, coupled problem on the spatial and temporal scales relevant to orogenic-gold systems is not currently possible within the timeframe of an exploration program. It is proposed that limiting the numerical models to the processes that represent the critical first-order controls associated with the system of interest can provide useful insight and predictive capabilities, so long as the chosen subset is appropriate to the system and scale under consideration. For orogenic-gold systems, the initial premise is that the critical controls on the localisation of ore-forming fluids, at the camp to deposit scale, are dilation, fluid flow and fluid focusing, and that these processes are controlled by rheological and permeability contrasts, the geometry of the system, and the stress regime at the time of mineralisation. Thus, numerical models of coupled deformation and fluid flow can be used to simulate orogenic-gold systems. However, it must be reiterated that these models do not consider thermal or chemical controls on mineralisation; hence the potential influence of these controls must be kept in mind when interpreting the results.

Two case studies are described which illustrate different approaches to ore-system modelling. The Stawell example demonstrates the use of a realistic 3D geometry to simulate a well-constrained conceptual geological model in the model-validation phase, and the subsequent application of the outcomes of these models to predictive exploration for similar deposits under cover. The Kundana example takes a more generic approach aimed at understanding why some fault intersection geometries in the Kundana field are mineralised, while others are not. A suite of models with varying geometries were run to identify those fault geometries with a higher probability of dilation and associated fluid focusing.

In both cases, a two-stage process is applied which involves reproducing the location of know mineralisation (the model validation phase), then applying the results to a predictive modelling phase aimed at ranking known targets or defining new prospective regions, targets or structures. Importantly, it is not necessary to prove a hypothesis beyond all doubt before applying it to exploration decision-making. The value proposition is based around identifying targets and detection methods with a higher probability of success than more empirical approaches. It is therefore critical that the results must be non-intuitive at some level. In both the Stawell and Kundana examples, it is argued that the target locations could not have been as well constrained on purely empirical grounds. In both cases, the higher-grade gold intersections would have been obtained eventually through closely spaced drilling, but at the expense of drilling many lower probability targets first.

Hydrothermal ore formation is the result of the interplay between many factors at many scales. Despite this, these examples demonstrate that a focus on first-order effects can deliver useful results. Clearly, in both cases, mineral deposition was also influenced by larger scale architecture, chemical gradients and the thermal regime. Nonetheless, useful conclusions were obtained without explicitly considering these aspects in the analysis. In other types of mineral systems or at different scales, chemical/temperature gradients and/or large-scale architecture might provide the first-order controls and should therefore be the focus of numerical modelling experiments for predictive targeting.

In conclusion, numerical modelling can be a useful experimental tool for the application of process thinking to mineral targeting in orogenic-gold systems, which can add significant value to the mineral exploration process.

Acknowledgements

We acknowledge the continued support of the pmd*CRC throughout this project. In particular, we thank Bob Haydon and Andy Barnicoat for their support and encouragement of the CRC's numerical modelling program. In addition, we acknowledge the advice, help and support of Kas DeLuca, Scott Halley, Greg Hall, Gerard Tripp, John Beeson, Marcus Willson, Susan Drieberg, Ian Neilson, Jon Dugdale, Chris Wilson, Tim Rawling, Jamie Robinson and John Miller, and the participation of Placer Dome (now Barrick Gold) and MPI Mines (now Perseverance Corporation Ltd.) in these projects. We would particularly like to thank one of the initial anonymous reviewers who took the time to provide detailed constructive criticism of the manuscript. This, along with the follow-up reviews, resulted in significant improvement to the paper.