Modelling of deformation around magmatic intrusions with application to gold-related structures in the Yilgarn Craton, Western Australia

Tectonophysics 526–529 (2012) 133–146

Contents lists available at SciVerse ScienceDirect

Tectonophysics journal homepage: www.elsevier.com/locate/tecto

Modelling of deformation around magmatic intrusions with application to gold-related structures in the Yilgarn Craton, Western Australia Y. Zhang a,⁎, A. Karrech a, b, P.M. Schaubs a, K. Regenauer-Lieb a, c, T. Poulet a, c, J.S. Cleverley a a b c

CSIRO Earth Sciences and Resource Engineering, PO Box 1130, Bentley, WA 6102, Australia School Civil and Resource Engineering, University of Western Australia, Crawley, WA 6009, Australia School of Earth and Geographical Sciences, University of Western Australia, Crawley, WA 6009, Australia

a r t i c l e

i n f o

Article history: Received 18 December 2010 Received in revised form 12 August 2011 Accepted 19 August 2011 Available online 31 August 2011 Keywords: Thermal continuum damage mechanics Rock deformation Intrusion Yilgarn Shear zone

a b s t r a c t This study simulates rock deformation around high temperature granite intrusions and explores how gold bearing shear zones near intrusions were developed in the Yilgarn, using a new continuum damage mechanics algorithm that considers the temperature and time dependent elastic-visco-plastic constitutive behaviour of crustal materials. The results demonstrate that strain rates have the most significant effects on structural patterns for both extensional and compressional cases. Smaller strain rates promote the formation of narrow high-strain shear zones and strong strain localisation along the flank or shoulder areas of the intrusion and cold granite dome. Wider diffuse shear zones are developed under higher strain rates due to strain hardening. The cooling of the intrusion to background temperatures occurred over a much shorter time interval when compared to the duration of deformation and shear zones development. Strong strain localisation near the intrusion and shear zone development in the crust occurred under both extensional and compressional conditions. There is always clear strain localisation around the shoulders of the intrusion and the flanks of the “cold” granitic dome in early deformation stages. In the models containing a pre-existing fault, strain localisation near the intrusion became asymmetric with much stronger localisation and the development of a damage zone at the shoulder adjacent to the reactivated fault. At higher deformation stages, the models produced a range of structural patterns including graben and half graben basin (extension), “pop-up” wedge structures (compression), tilted fault blocks and switch of shear movement from reverse to normal on shear zones. The model explains in part why a number of gold deposits (e.g. Wallaby and Paddington deposits) in the Yilgarn were formed near the flank of granitecored domes and deep “tapping” faults, and shows that the new modelling approach is capable of realistically simulating high strain localisation and shear zone development. Crown Copyright © 2011 Published by Elsevier B.V. All rights reserved.

1. Introduction In recent years, numerical modelling has emerged as a powerful supplementary tool for studying mineralisation processes and assisting mineral exploration (e.g. Garven et al., 2001; Ord et al., 2002; Potma et al., 2008; Schaubs et al., 2006; Sorjonen-Ward et al., 2002). This is due to its computational capability to integrate data and efficiently simulate the interactions between rock deformation, fluid flow, thermal transport and chemical reactions during mineralisation, the four key processes generally leading to ore formation (e.g. Hobbs et al., 2000; Zhang et al., 2003; Zhao et al., 1998). The late Archaean Yilgarn Craton of Western Australia, one of the most important terrains hosting orogenic gold deposits in the world (e.g. Cassidy et al., 1998; Goldfarb et al., 2001; Groves et al., 2000), represents an excellent terrain for the application of numerical modelling

⁎ Corresponding author. Tel.: + 61 8 6436 8626; fax: + 61 8 6436 8555. E-mail address: [email protected] (Y. Zhang).

and predictive targeting. The main reason is that gold mineralisation in the Yilgarn is characterized by strong structural controls, extensive involvement of hydrothermal fluids and the timing of mineralisation being coeval with metamorphism and closely linked to intrusive events (e.g. Blewett et al., 2010a,b; Cassidy and Hagemann, 2001; Davis et al., 2010; Groves et al., 2000; Weinberg and van der Borgh, 2008). Numerical modelling has already been used in the Yilgarn to study mineralisation processes and assist mineral exploration (e.g. Holyland and Ojala, 1997; Potma et al., 2008; Sorjonen-Ward et al., 2002). The majority of the previous coupled deformation and fluid flow models for the Yilgarn and also those for other regions in the world (e.g. Ord et al., 2002; Schaubs et al., 2006; Zhang et al., 2006, 2007b) are based on elastic-plastic constitutive behaviours (strain-rate and temperature independent) and focused on the deformation effects of pre-existing structures (mostly faults) on fluid flow. The initiation and development of structures under thermal influence and rate dependence have not been included in such upper-crust scale studies with a mineralisation focus. Rate-dependent deformation behaviours and thermal effects could be very important for the realistic modelling of rock deformation and structural evolution, particularly for the Yilgarn. Recent structural and

0040-1951/$ – see front matter. Crown Copyright © 2011 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2011.08.013

134

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

seismic studies (e.g. Blewett et al., 2010a,b; Davis et al., 2010; Goscombe et al., 2007; Henson et al., 2007) show that gold mineralisation and the development of Au-bearing shear zones may be structurally and temporally linked to magmatic intrusions in the Yilgarn. These intrusions include mafic granite intrusions around 2.665 Ga, syenites around 2.665–2.650 Ga or more extensive low-Ca granite intrusions at/from 2.65 Ga. Several world-class gold deposits and gold-hosting shear zones are observed to have occurred above or on the margins of intrusions (granite domes) (Blewett et al., 2010b; Czarnota et al., 2010a; Henson et al., 2007, 2010); the contacts between granite intrusions/domes and greenstones are often interpreted as faults and shear zones (Blewett et al., 2010a). To simulate the effects of high temperature granite intrusions on the development of high strain location (e.g. shear zones), the inclusion of rate-dependent behaviours and thermal contribution in a numerical model is important. Regenauer-Lieb et al. (2006) developed a fully-coupled thermal-mechanical model and simulated the reduction of lithospheric strength as a consequence of the feedback between deformation, shear-heating and temperature dependence of flow laws. It was shown that such strength reduction can result in the development of mid-crust detachment faults. This approach has been further advanced recently by Karrech et al. (2011a) in the development of a continuum damage mechanics algorithm that considers the elastic-visco-plastic constitutive behaviour of crustal materials experiencing continuum damage with dependence on time, temperature, pressure and water content. In this study, we apply this new continuum-damage-mechanics algorithm to a generic architectural model relevant to the Yilgarn, which contains a high-temperature magmatic intrusion. Our aims are two fold: (1) explore the applicability of the new theories and algorithm of Karrech et al. (2011a) to a key mineralisation terrain in terms of simulating strain localisation; (2) gain some understanding of rock deformation behaviour around high temperature intrusions (domes) in the Yilgarn and hence provide a generic reason why several world-class gold deposits occurred above or on the margins of intrusions (e.g. Blewett et al., 2010b; Henson et al., 2007). This modelling study does not consider the effects of fluid flow or chemical reactions, and as such is unable to directly predict gold mineralisation. 2. Modelling methodology The numerical model developed herein uses the damageable constitutive behaviour of geomaterials as described by Karrech et al. (2011a). In this section, we provide a brief description of the mathematical formulation and emphasize the key aspects of the underlying this continuum damage model. The constitutive model is based on the nonlinear deformation of materials, which takes into consideration a coherent scheme of energy exchange in terms of storage and dissipation. This scheme relies on the expressions of the Helmholtz free energy and Clausius–Duhem equations, respectively. By applying the principle of maximum dissipation, this scheme results in constitutive relations and flow rules which govern the weakening processes (Karrech et al., 2011b; Regenauer-Lieb et al., 2010). 2.1. Energy storage and dissipation The processes which govern material behaviours are generally controlled through state variables. The relevant state variables considered in this paper are deformation, temperature and damage. They respectively measure the kinematics, thermal transfers and failure of the considered bodies. The material kinematics is governed by a deformation gradient of the form F = ∇ ϕ and ϕ is a mapping function from the initial configuration to the current one. This results in heat exchanges which are driven by temperature, T. These processes trigger materials degradation which is described by a scalar, D, and a set of other dissipative parameters which are summarized by ξ. These quantities can be used to describe the Helmholtz free energy of the system as follows: ψ(F, D, T, ξ). Applying the first principle of

conservation of momentum and energy results in the following equations (Regenauer-Lieb et al., 2010):


divðσÞ þ ρf ¼ 0 in ρcpT˙ ¼ −divðqÞ þ σ : d þ r

ð1Þ

Where, σ is the Cauchy Stress, ρ is density, f is a body force, cp is a heat capacity, q is a heat flux and r is a heat source. d in denote the dissipative part of the deformation rate whichcan be expressed  in terms ˙ −1 þ F−t F˙ t =2. Solving of gradient of deformation as follows:d ¼ FF Eq. (1) requires the Dirichlet and/or Neumann boundary conditions as well as constitutive relations. Applying the second principles of thermodynamics and using the relationships between specific entropy, free energy, and specific internal energy, the Clausius–Duhem inequality is expressed as follows: σ : d−

dψ q − :gradT ≥0 dt T

ð2Þ

The heat flux is expressed as q = − kT gradT, where kT is a thermal conductivity. Independency of the different processes along with the definitions of Helmholtz free energy and the Clausius–Duhem inequality result in the following relationships: σ¼

dψ dψ ; and Y ¼ : de dD

ð3Þ

Where e is an appropriate strain measure (Karrech et al., 2011c). In case of linear elasticity, the intrinsic (non thermal) thermodynamic forces acting the material obey to: 8 > > ˙ > <σ ¼

   νð1−DÞE ð1−DÞE  in in tr d−d Ι þ d−d þ Ωσ−σΩ ð1 þ νÞð1−2νÞ ð1 þ νÞ " !2 # : ð4Þ 2 σ 2 σH > eq > > ð 1 þ ν Þ þ 3 ð 1−2ν Þ Y ¼ : 2 3 σeq 2Eð1−DÞ The first term in the above equation represents the co rotational rate of Cauchy stress and the second denotes the thermodynamic force of damage. The constants E and v represent respectively the Young modulus and Poisson's pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiratio of the undamaged material. σH = tr(σ)/3 and σeq ¼ 3s : s=2 are stress invariants and s is the deviatoric stress. The dissipative part of the stretching rate is calculated using a flow rule of the form: in

d ¼λ

∂g : ∂σ

ð5Þ

Where g is a Drucker–Prager plastic potential: g(σ) = σeq − ασH. A similar form was used to describe the limit of elasticity as follows: f(σ) = σeq − βσH − κ − φ [− 1]. In the above equations the parameters α and β are directly related to the friction and dilation angles. The viscous effects are considered through the function φ [− 1] also known as “overstress”. Its inverse is expressed in an additive decomposition of dislocation and diffusion creep terms as follows:     Q σ Q n φ ¼ ε˙ dis þ ε˙ dif ¼ A1 σ exp − 1 þ A2 exp − 2 : ð6Þ γ RT RT Ai and Qi are the pre-coefficients and activation energies of the two mechanisms, respectively. In order to describe the damage evolution we proceed in a similar manner as for the inelastic process. We first define a damage potential: gY =((1−D)− {n + 1} −1)Y+ℵ(y), where the first term refers to the damage evolution and the second describes the damage nucleation. This expression is then used to calculate a flow rule of the form: ∂g D˙ ¼ Y : ∂Y

ð7Þ

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

The problem under consideration consists in solving Eq. (1) using the constitutive relations (Eqs. (2) and (4)) while taking into account the flow rules (Eqs. (3) and (5)). This problem is solved numerically using Abaqus as well as a fortran subroutine routine which integrates Eqs. (4)–(7).

135

Effects of temperature are included through two different coupling mechanisms. The first is thermal expansion which accounts for the deformation due to variation of temperature and the second is shear heating which accounts for the change of temperature due to dissipation. Given these coupling effects, the drop of temperature in any point of the model will automatically generate a proportional material shrinkage.

2.2. Strain localisation due to damage 3. Modelling application to the Yilgarn A novel element in the continuum damage mechanics approach introduced by Karrech et al. (2011a) is that it extends the damage formulation to include materials failure when viscous effects are active. This is generally the case in the brittle–ductile transition zone and below. The Fig. 7 of Karrech et al. (2011a), illustrates the behaviour of materials in this transition zone under various loading velocities. In accordance to Ashby et al. (1979), three different fracture mechanisms can take place: (i) brittle fracture: characterized by elastic behaviour without significant dissipation prior to failure, (ii) creep fracture: which requires rate dependent dissipation before failure, and (iii) ductile fracture: involving significant rate independent dissipation before failure. Slow deformation can take place entirely in the elastic regime and it can have no significant dissipation prior to failure and may consequently be classified as a brittle fracture. Since our formulation includes all three rheological bodies, elastic, viscous and plastic, instability caused by damage can occur in all three regimes and all three fracture mechanisms are possible. The particular responses are attributed entirely to the coupled elasto-visco-plastic behaviour embedded in the framework. These processes are time dependent, that is, failure does not depend only on materials strength but also on the rates of loading. Therefore, the resultant geometrical features such as width of shear zones and their orientations are controlled not only by the directions of stresses but also by the rates of loading as well as the distributions of temperatures. A novel aspect of our formulation is that we consider the importance of creep fracture. Creep fracture operates very efficiently at low loading rates but much less efficiently at high loading rates. This is exhibited as rate-dependent strain hardening behaviour with the increase of loading or strain rates as described by Karrech et al. (2011a). A precursor to this surprising mechanism has recently been discovered to be important in ductile shear zones in granites (Fusseis et al., 2009). The precursor to creep fracture is supported by a process called creep cavitation, i.e. the generation of micro-pores on grain boundaries linked to slow creep deformation. Its existence has been postulated earlier based on circumstantial evidence of volatile transfer through the ductile part of the lithosphere (Regenauer-Lieb, 1999). Direct geological evidence for ductile fracture in the presence of melts has been supplied by (Weinberg and Regenauer-Lieb, 2010). The addition of these two new mechanisms, creep and ductile fracture to a forward model of crustal deformation obviously leads to much more efficient localisation. However, it may require further geological evidence before being widely accepted since it directly challenges geological intuition that high strain rate favours brittle, localised deformation while low strain rates facilitate distributed viscous deformation. Evidence for the importance of localisation by slow processes such as creep and ductile fracture rely for the time being mainly on forward prediction of the model formulated by Karrech et al. (2011a). It predicts a quadratic energy with respect to deformation, before failure takes place. However, unlike non damageable materials which exhibit a plateau like variation of energy in the plastic regime, the approach shows how stored energy decreases with respect to loading. This response, which is directly related to degradation, explains how materials dissipation due to damage is fundamentally different to materials dissipation due to dislocations or diffusion. The models also reveal that visco-plastic dissipation is higher in non-damaged structures. Therefore the necessary cost of deforming geological materials with classical approaches is highly overestimated.

3.1. Architecture of generic models for the Yilgarn The Yilgarn Craton has a long, complex geological history involving at least 6 deformation events and multiple magmatic and metamorphic episodes (e.g. Blewett et al., 2010a; Czarnota et al., 2010b; Goscombe et al., 2007), which led to the development of a regional architecture. The geology of the region is characterized by an extensive greenstone succession overlying a granite domain, together with younger sedimentary rocks confined to fault-bounded late basins. The greenstone (tholeiitic and komatiitic mafic-ultramafic rocks plus calcalkaline felsic volcanics and clastic sedimentary rocks) is the product of a series of volcanic episodes which occurred between ~2710 and ~2685 Ma, approximately associated with the M1 metamorphic event (Goscombe et al., 2007; also see Blewett et al., 2010a). The bulk (60%) of the granitic domain is represented by high-Ca granites emplaced approximately between 2685 and 2665 Ma, coinciding with the M2 regional metamorphic event (Goscombe et al., 2007). Later mafic granites (between ~2665 and 2650 Ma), low-Ca granites (between ~2650 and ~2620 Ma) and syenites (~2665–2650) were intruded into the high-Ca granite and greenstone rocks. These later intrusions are generally associated with major gold mineralisation phases in the region (e.g. Blewett et al., 2010a). The present generic models are designed to simulate the structural framework for a particular point of time during the geological history for the Yilgarn, that is, following the formation of the greenstone and high-Ca granite successions and immediately after the emplacement of a hot magmatic intrusion (the intrusion can be related to mafic granite, low-Ca granite or syenite intrusions). The key questions here include: 1. How do the high temperature intrusions affect the deformation of host rocks and, 2. what are the key factors controlling the development of shear zones? As such, a simple 2D architecture from a cross section view has been adopted for the models (Fig. 1). The section is 40 km long and 10 km deep, containing a greenstone unit, a granitic unit and a hot granite intrusion. The granite unit has a dome-shaped geometry in the central part of the model, reflecting the shape of structural features of the felsic upper crust (granitic domain) seen in the Yilgarn (e.g. Henson and Blewett, 2006). The hot granite intrusion is located within a domal crest of the cold granite unit. Considering the proposal of Blewett et al. (2010a) that fault structures had most likely already developed before later intrusions, the effects of a pre-existing fault has been explored in some models (Fig. 1). From a generic modelling point of view, it is necessary to explore a range of geometry scenarios for systematic variations of geometries from simple to complex (e.g. no fault versus with faults). 3.2. Model properties The elastic-plastic, thermal and viscous creep properties are given in Table 1. The elastic-plastic and thermal properties are based on the average parameters for mafic rocks and granite in the literature (Turcotte and Schubert, 1982), with specific heat parameters from an engineering data set (http://www.engineeringtoolbox.com/specific-heat-solidsd_154.html). The parameters for dislocation creep are based on the

136

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

Fig. 1. (a) Geological architecture of the numerical models. (b) Numerical mesh for a portion of the model (see dash-line outlined area in a). (c) Initial temperatures of the model. The intrusive body has the highest temperature of 973.17 K (700 °C).

results of Ji et al. (2003) for diabase and granite. A set of homogeneous diffusion creep parameters (based on data from Ranalli, 1997, 2000) are adopted in the model to simulate the homogeneous contribution of diffusion creep from all the rock units. This set-up simulates the domination of dislocation creep together with elastic-plastic response in the upper crust levels and emphasizes the different dislocation creep behaviours of different rock types. In the models containing a pre-existing fault, the fault is simulated as a narrow, soft elastic-plastic zone with Young's modulus = 3.5 × 1010 Pa, Poisson's ratio = 0.25 and plastic strength = 1.5 × 107 Pa for simplicity. 3.3. Boundary conditions A recent structural study (Blewett et al., 2010a) showed that the Yilgarn Craton experienced subsequent deformation events involving extension and contraction or transpression (i.e. D3 onward) after the formation of greenstone and high-Ca granite units, and that these correlate with major gold events in the region. It is therefore necessary to explore both extensional and contractional deformation scenarios in the present models. Divergent and convergent velocity conditions are applied to the two vertical edges of the model (Fig. 1), >respectively, to create extensional and compressional (shortening) deformation cases. The bottom boundary of the model is not allowed to move in the vertical

direction but can slip in the horizontal direction. The values of the velocities are determined based on the required strain rate. A model of present-day plate motions and plate boundary deformation (Kreemer et al., 2003), which included 3000 geodetic velocities, suggested that strain rates for active plate margins are in the range of 1.0 × 10−15 to 1.0 × 10−13 s −1. On the basis of the predictions above, we have adopted the velocity boundary conditions for several different strain rate cases: 1) 0.5×10−15 s−1; 2) 1.0×10−15 s−1; 3) 0.5×10−14 s−1; and 4) 0.5×10−13 s−1. Note that the deformation velocities at the two vertical edges of the model are fixed at the initial values and hence the strain rates listed above are bulk or average strain rates. Thermal boundary conditions for the models, include: 1) a fixed surface temperature of 298.18 K (25 °C) for the top of the model; 2) a geothermal gradient of 30 K/km, leading to a temperature of 598.18 K (325 °C) for the model base; 3) an intrusive temperature of 973.18 K (700 °C) for magmatic intrusion – this intrusive temperature is allowed to cool with time; and 4) a thermal flux of 35 mW m −2 is applied to and maintained at the model base throughout simulation. The application of the homogeneous thermal flux above along the base of the model is a simplified approach. To maintain a geothermal gradient of 30 K/km throughout simulation, a carefully-calibrated and varying thermal flux (along the model base and with simulation time) would have to be implemented due to the presence of in-homogenous

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

137

Table 1 Key material properties of the models. Model parameters Elastic-plasticity Density (kg.m−3) Young's modulus (Pa) Poisson's ratio Plastic strength (Pa) Thermal parameters Thermal conductivity (W m−1 K−1) Specific heat (J kg−1 K−1) Thermal expansion coefficient (K−1) Dislocation creep parameters Universal gas constant (J mol−1 K−1) Activation energy (J mol−1 K−1) Prefactor of power-law equation (MPa−n s−1) Stress component (n) Diffusion creep parameters Activation energy (J mol−1 K−1) Prefactor of power-law equation (MPa−n s−1) Stress component (n)

Greenstone

Granite

Intrusion (hot)

Fault

2950 7.0 × 1010 0.25 4.0 × 107

2650 6.0 × 1010 0.25 2.0 × 107

2600 1.0 × 1010 0.25 0.5 × 107

2650 3.5 × 1010 0.25 1.5 × 107

2.5 840 2.4 × 10−5

3.0 790 2.4 × 10−5

3.0 1000 2.4 × 10−5

2.8 1200 2.4 × 10−5

8.3144 3.519 × 105 6.12 × 10−5 4.97

8.3144 2.06 × 105 5.74 × 10−7 3.9

8.3144 2.06 × 105 1.0 × 10−3 3.9

1.5 × 105 4.8 × 10−4 1.1

1.5 × 105 4.8 × 10−4 1.1

1.5 × 105 4.8 × 10−4 1.1

model internal structures and development of topography at model top as the result of bulk deformation.

4. Model results 4.1. Extensional case at a strain rate of 1.0 × 10 −15 s −1 Our first numerical simulation investigates an extensional scenario with a strain rate of 1.0 × 10−15 s −1, deformed to approximately 10% bulk extension. Fig. 2 presents the strain patterns for four stages during the simulated deformation history. In the very early deformation stage (0.4% bulk extension or 0.17 million years after intrusion), strain localisation mainly occurred at the two shoulder locations of the hot intrusive body (Fig. 2). There is also less prominent strain localisation along the boundary of the intrusion and some localisation at the flanks of the “cold” granite dome (see Fig. 1a for model geometries). These locally-higher strains predominantly reflect plastic-viscous strains, whereas the most parts of the model (very low strain) are still in the elastic mode. By the 0.69% bulk extension stage (1.0 million years after intrusion), strain continued to be localised at the shoulder locations and boundary of the intrusion (see strain magnitudes in Fig. 2b), but strain localisation along the two flanks of the dome was more intensive, with high strain values and clear shear zone formation along these flanks. With continued strain localisation at 2.11% bulk extension (1.83 million years after intrusion), structural patterning (Fig. 2c) became more obvious and is characterised by: (1) two shear zones along the dome flanks continued to grow and intersected in the greenstone unit above the granite dome; (2) two new shear zones were formed from earlier strain localisation at/along the intrusion shoulders and boundary, extending into the greenstone unit; (3) several new smaller scale shear zones were developed within the intrusion. By 10% bulk extension (4.03 million years after intrusion), the structural patterns above remained in place but were enhanced (Fig. 2d). Two pairs of shear zones exhibit large displacement, and this final geometry of the model suggests the formation of well-developed extensional double-graben basins bounded by newly developed shear zones. For comparison, a reference model under the same boundary conditions but without the intrusion (see Fig. 1a) was performed. The final strain patterns at ~ 10% bulk extension (Fig. 2e) shows the spacing between the developed shear zones is wider than in the case with the intrusion (Fig. 2d), due to the absence of the intrusion and its influence. Examination of strain evolution in this model indicates that the central pair of shear zones was initiated entirely along the flanks and apex of the dome. They then propagated upward and downward

with increased bulk deformation and became straightened at the end of ~10% bulk extension (hence the highest strain zones shifted away from the edge of the dome). Plots of temperature variations through several stages of the simulation with the intrusion (Fig. 3) reveal quite rapid cooling of the intrusion. At 0.4% bulk extension or 0.17 million years after intrusion, the temperature of the intrusive body has already cooled by as much as ~ 365 °C, with remnant high temperatures (anomalous higher temperatures relative to the background geothermal state) around the bottom boundary of the intrusion. The intrusive temperatures were almost cooled to the background geothermal state by 0.69% bulk extension or 1.0 million years after intrusion emplacement (Fig. 3b). Later stage temperature variations (Fig. 3c and d) are predominantly due to the change of crustal thickness as a result of extensional deformation.

4.2. Compressional case at a strain rate of 1.0 × 10 −15 s −1 For the compressional case, the model produced a different type of geological architecture (Fig. 4), but with some similar structural characteristics to the corresponding extensional model (Figs. 2 and 3). Early strain localisation (0.6% bulk shortening and 1.0 million years; Fig. 4a) still occurred at the shoulder areas of the intrusion (also the top boundary of the intrusion) and the flanks of the “cold” granite dome. There is a clear link (via relatively higher strain zones) between the maximum strain spots at the intrusion shoulders and along the flanks of the dome. At 1.03% bulk shortening (1.31 million years after intrusion), two obvious shear zones developed, extending into and intersecting the greenstone unit. These shear zones display reverse movement, resulting in some thickening of the crust. From the 3.0% bulk shortening stage (Fig. 4c) and finally to the 10% bulk shortening stage (Fig. 4d), the two shear zones along the flanks of the dome continue to be enhanced, and two new shear zones are developed more distally. The latter shear zones seem to be predominantly evolved from the early maximum strain localisation at the shoulder areas of the intrusion plus strain localisation along the bottom boundary of the intrusion. All these shear zones exhibit a large amount of reverse movement and led to the large thickening of the crust and the development of a wedge-shaped “pop-up” or “flower”-type structure. This is in clear contrast to the results of the extensional model, where normal shearing movement and crustal thinning occurred. Temperature distributions again show the cooling of the intrusion is rather rapid (Fig. 5). At the 0.6% bulk shortening (1.0 million years after intrusion) stage (Fig. 5a), the temperatures mainly reflect the background

138

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

Fig. 2. Development of shear strain patterns for extensional deformation models at an overall strain rate of 1.0 × 10−15 s−1. a, b, c and d are for the bulk extension stages of ~ 0.4%, ~ 0.69%, ~ 2.11% and ~ 10.0%, respectively. e is the final pattern of strain at ~ 10% bulk extension for a reference model that only contains the greenstone and “cold” granite units but without the intrusion (see Fig. 1a).

geothermal state with probably only some remnant effects from the intrusion. A reference model without the intrusion under compressional conditions was also performed and the strain pattern at ~10% bulk shortening is given in Fig. 4e. Similar to the extensional case (Fig. 2), the major effect of the exclusion of the intrusion in the model is reflected by the change in shear zone spacing. The development of the central pair of shear zones is entirely controlled by the flanks of the dome and shows a relatively smaller spacing, in comparison with the model containing the intrusion (Fig. 4d). Two outer shear zones developed further away from the central shear zones (spacing increased). These changes are due to the absence of the intrusion and its effects on the initiation and development of shear zones in the model.

4.3. Effects of strain rate variations 4.3.1. Extensional case To understand the effects of strain rates on deformation patterns under extensional conditions, three more extensional models with different strain rates (0.5× 10−13, 0.5 × 10 −14 and 0.5 × 10 −15 s−1) are investigated and the results are compared with that for the strain rate of 1.0 × 10−15 s −1 below (Fig. 6a–d). A key result from this set of models is that shear zones become more diffuse with increasing strain rate, in company with the development of fewer shear zones. Increasing the strain rate from 1.0× 10−15 s −1 (Fig. 6c) to 0.5 × 10−14 s −1 (Fig. 6b) led to quite different strain distribution patterns and shear zone geometries. Under the higher strain rate, shear

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

139

Fig. 3. Temperature distributions in the extensional model with a bulk strain rate of 1.0 × 10−15 s−1 shown in Fig. 2. a, b, c and d are for the bulk extension stages of ~ 0.4%, ~ 0.69%, ~ 2.11% and ~ 10.0%, respectively.

zones became more diffuse with lower strain values (Fig. 6b) relative to the lower strain-rate case (Fig. 6c). These more diffuse shear zones are predominantly initiated from high strain localisation along the flanks of the high temperature intrusive body during the early stage of the deformation history. There is no major strain localisation along the flanks of the “cold” granite dome and hence, the central pair of shear zones which developed in the lower strain rate case (Fig. 6c) is missing here (Fig. 6b). In the case where strain rate was further increased to 0.5 × 10−13 s −1 (Fig. 6a), strain in the “cold” greenstone and granitic units is widespread (more homogeneously distributed), lacking clear shear zones. Strain localisation is almost entirely confined to the hot intrusion, exhibited as small shear zones along the boundary of or within the intrusion. Decreasing the strain rate to 0.5 × 10−15 s −1 resulted in more profound changes in the strain pattern (Fig. 6d). Shear zones in this model became remarkably asymmetric. More importantly, strain is more localised, leading to more shear zones developed here than in the 1.0 × 10 −15 s−1 case (compare Fig. 6d with Fig. 6c; also see Fig. 2). Although there is early strain localisation along the flanks of the intrusion and granitic dome (similar to the patterns shown in Fig. 2b), the first pair of emerging shear zones (F1a and F1b in Fig. 6d) are partially along the flanks of the “cold” granite dome.

With the continuation of bulk extension, F1a became more developed and accumulated greater strain than F1b. This led to the initiation and development of F2 due to the need to accommodate the normal movement along F1a. Further extension resulted in the initiation and development of three more shear zones, F3a, F3b and F3c, while earlier shear zones continued to accumulate more shear strain. The final structural outcome at 10% bulk extension is a complex shear zone network confining tilted “fault” blocks and asymmetric graben basins. 4.3.2. Compressional case Three more compressional models with different strain rates (0.5 × 10 −13, 0.5 × 10−14 and 0.5 × 10−15 s−1) were also conducted. The comparison of the modelling results between these models and the model with the strain rate of 1.0 × 10−15 s−1 also clearly illustrates the effects of strain rates on the patterns of strain distribution and shear zone geometries (Fig. 6e–h). Similar to the extensional cases, the modelling results show that increasing strain rate also led to increasingly diffuse shear zones. When the strain rate is increased from 1.0×10−15 s−1 to 0.5×10−14 s−1 (Fig. 6f), shear zones became clearly more diffuse. A wedge-shaped “pop-up” structure and inhomogeneous crustal

140

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

Fig. 4. Development of shear strain patterns for compressional models at a bulk strain rate of 1.0 × 10−15 s−1. a, b, c and d are for the bulk shortening stages of ~ 0.6%, ~ 1.03%, ~ 3.0% and ~ 10.0%, respectively. e is the final pattern of strain at ~ 10% bulk shortening for a reference model that only contains the greenstone and “cold” granite units but without the intrusion (see Fig. 1a).

thickening is still developed, but the main reverse shear zones bounding the “pop-up” structure are significantly wider than the corresponding shear zones in the model with a lower strain rate of 1.0 × 10−15 s−1

(Fig. 6g; also see Fig. 4). The initiation of these diffuse shear zones are mainly governed by early strain localisation in and around the intrusion (i.e. shoulder and boundary sites). There is strain localisation along

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

141

Fig. 5. Temperature distributions in a compressional model at a bulk strain rate of 1.0×10−15 s−1. a, b, c and d are for the bulk shortening stages of ~0.6%, ~1.03%, ~3.0% and ~10.0%, respectively.

the flanks of the dome in the early deformation stages, but these high strain zones were overprinted by the diffuse shear zones and hence did not develop into individual shear zones. When the strain rate is further increased to 0.5 × 10 −13 s−1 (Fig. 6e), strain localisation in the “cold” rocks became even more diffuse so that strain in these rock units is near homogeneous and there is a lack of distinct shear zones. Major strain localisation only occurred within the hot intrusion and this is expressed as the development of small shear zones along the intrusion boundary and within the intrusion. Decreasing the strain rate to 0.5 × 10 −15 s −1 resulted in the increased shear localisation (Fig. 6h). Shear zone patterns became asymmetric, similar to the extensional case (Fig. 6d). Examination of shear zone initiation and development sequence reveals that shear zones F1a and F1b were first initiated in a manner similar to the patterns illustrated in Fig. 2a and b. When reverse movement or thrusting along F1a became much greater than F1b; propagation of the latter essentially stopped. The thrusting of F1a led to the formation

of F2 as a “back-thrust” accommodate the bulk strain. This led to the development of an asymmetric, wedge-shaped “pop-up” structure (Fig. 6h). The formation of shear zones F3a and F3b occurred later in the deformation history. 4.4. Effects of a pre-existing fault Two models with a pre-existing fault were performed for extensional and compression deformation conditions, respectively. The fault in the model is simulated as a narrow, soft elastic-plastic fault zone (i.e., not as a discrete fault plane) which is a simplification of the complex features of natural faults. Small faults in nature can occur as discrete planes, but it has been increasingly recognized that faults (particularly large faults) often occurred as fault zones (see review by Wibberley et al., 2008) which are expressed as high strain or damage zones. In the literature, geologists sometimes call such high strain zones equally as either faults or shear zones (e.g. Wang et al., 2007).

142

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

Fig. 6. Plots of shear strain patterns at 10% bulk extension or shortening for eight numerical models at four different strain rates, illustrating the effects of strain rates. a, b, c and d are for extensional models at strain rates of 0.5 × 10−13, 0.5 × 10−14, 1.0 × 10−15 and 0.5 × 10−15 s−1, respectively. e, f, g and h are for compressional models at strain rates of 0.5 × 10−13, 0.5 × 10−14, 1.0 × 10−15 and 0.5 × 10−15 s−1, respectively. F1 to F3 in d and h identify major shear zones developed in the models (numbers represent the order that shear zones emerged, i.e. 1 is the earliest).

It is noted that in the early deformation stages under both conditions, strain localisation around the intrusion is clearly asymmetric (Fig. 7a and d), in contrast to the results of the models without a pre-existing fault. Greater strain localisation occurred near the left shoulder area of the intrusion, which is adjacent to the pre-existing fault, than the right shoulder area. This is exhibited as the development of a higher strain zone linking the intrusion to the fault in very early deformation stage. The differences of shear zone development patterns between the extensional and compression cases are described below. 4.4.1. Extensional case Early during extensional deformation (0.66% bulk extension), the F1 normal shear zone first developed in a manner similar to the antithetic faulting mechanism, along with large normal movement along

the pre-existing fault (Fig. 7a). There are also higher strains localised near the left shoulder of the intrusion. By ~4% bulk extension (Fig. 7b), normal movement along the pre-existing fault increased, leading to the further development of F1. Initially higher strains at the left shoulder of the intrusion did not develop into clear shear zones by this stage. The development of the structures above led to the formation of asymmetric wedge-shaped, tilted fault block and half-graben geometries. It is also noted that the F2a and F2b shear zones were initiated next, but these shear zones remain as more diffuse (i.e. lower strains than other more dominant shear zones such as F1) by the 10% bulk extension stage (Fig. 7c). The F3a, F3b and F3c shear zones developed with the continuation of bulk extension; where F3c is the result of the enhancement and propagation of earlier higher strains near the left shoulder of the intrusion in later deformation

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

143

Fig. 7. Development of shear strain patterns in numerical models containing a pre-existing fault. a, b, and c are for an extensional model at 0.66%, 3.96% and 10.0% bulk extension, respectively. d, e, and f are for a compressional model at 0.6%, 3.9% and 10.0% bulk shortening, respectively. The strain rate for the both models is 1.0 × 10−15 s−1. F1 to F3 in a, b and c, and F1 to F4 in d. e and f identify major shear zones developed in the models (numbers represent the order that shear zones emerged, i.e. 1 is the earliest).

stage. The final architecture of the model at 10% bulk extension includes an extended and thinned crust characterized by the formation of multiple normal shear zones/faults and an asymmetric multiple half graben basin (Fig. 7c). This model shows that when a pre-existing fault is included, strain localisation and structural development are dominantly controlled by movement along the fault, and the boundaries of the dome and intrusion only play a minor role in localising strains.

regime is still compressional or convergent. Such a switch might be partially attributed to the result of the exaggerated topography developed in the absence of erosion in the model. F5a and F5b were the last to develop in the sequence, and in particular, the initiation of F5a is the result of the thrust loading from the “pop-up” block between the initial fault and F2 (Fig. 7f).

5. Discussion 4.4.2. Compressional case Strong strain localisation near the left shoulder of the intrusion led to the development of a well defined shear zone, F1, at the 0.6% bulk shortening stage (Fig. 7d), providing an early damage link between the intrusion and the initial fault. The F2 shear zone then developed as a “back-thrust” shear zone (Fig. 7e), in the similar position as the main normal shear zone in Fig. 7b, due to reverse faulting or thrusting along the initial fault. This led to the development of an asymmetric wedge-shaped “pop-up” block, with greater elevation at the initial fault side. F3 then developed near the right corner of the “pop-up” block as the result of continued reverse faulting along the initial fault, and a secondary “pop-up” wedge block was formed between the initial fault and F3 (Fig. 7e). F4a and F4b were initiated next (Fig. 7f), while F1 continued to accumulate strain and grow (F4b and F1 are clearly separate shear zones at the time of initiation although they are almost amalgamated at the final deformation stage). It is interesting to note that the movement of F3 and F4b switched to normal or extensional shearing, even though the overall deformation

5.1. Strain rate, intrusion cooling rate and structural style The models presented here involve several different strain rates ranging from 0.5 × 10 −15 to 0.5 × 10 −13 s −1, but these represent average strain rates for a simulation. Local strain rates within high strain zones will be much higher than the rates for less deformed areas outside of the shear zones in the same simulation. The modelling results show that varying the strain rate has the most significant effect on the development of broad structural patterns for both extensional and compressional cases (Figs. 6 and 7). Smaller strain rates (0.5 × 10 −15 and 1.0 × 10 −15 s −1) promote the development of well-defined shear zones with highly localised strain (Fig. 6c, d, g and h), together with strong strain localisation along the flanks or shoulder areas of the intrusion and cold granite dome (See Figs. 2 and 4). In the lowest strain rate models (Fig. 6d and h), asymmetric structural patterns emerge, which reflect only a weak influence from the dome and intrusion. These are in clear contrast with the structural patterns of

144

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

diffuse shear zones (wider and lower strains) in the models with higher strain rates of 0.5 × 10 −14 and 0.5 × 10 −13 s −1 (Figs. 6a, b, and 7e, f). In particular, strain localisation in the highest strain rate models (Fig. 6a and e) occurred predominantly within the hot intrusion with near homogenous deformation in the upper crust (greenstone and granitic units). The changes of modelled strain localisation patterns with changes in strain rates are due to the strain hardening behaviour of the coupled elasto-visco-plastic rheology of the models at higher strain rates, which is associated with the change of the efficiency of the creep fracture mechanism with the change of strain rate, as discussed in Section 2 above and fully described in Karrech et al. (2011a). As strain rates increase, such materials display progressive hardening and hence more diffuse shear zones develop. At this point it is important to emphasise the new elements of our models described in Section 2.1, in particular the roles of stored and dissipated energy for damage generation. Our models allow a homogenized version of creep and ductile fracture modes, which are extremely efficient at low strain rates. If the strain-rates are increased the two mechanisms that rely on damage accrued by dissipation are effectively suppressed. The inability of strain weakening by creep and ductile damage, for high loading rates imply that our model produce a somewhat counterintuitive material behaviour with strain hardening at such rates. The analysis of present-day geodetic data (Kreemer et al., 2003) suggests that the strain rates of ~10−15 s−1 represent the upper bound of intra-plate strain rates or the lower bound of active plate margin strain rates. Under such strain rates, a 20% bulk extension or shortening takes about 6 million years to complete. This time scale is well within the durations of several prominent orogenic or tectonic events for the earth. For the Yilgarn, at least 6 deformation events are involved, and each of these events occurred over million years (e.g. ~10 million years for D3; see Blewett et al., 2010a). Other examples for some other regions in the world are: 1) the Mount Eclipse Movement of the Alice Springs Orogeny, Central Australia (~50 million years from ~350 to ~298 Ma; Dunster et al., 2007); 2) the India-Eurasia plate collision zone (16 million years from ~60 to ~44 Ma for “soft” collision and 40+ million years from ~44 Ma to present for “hard” collision; see Lee and Lawver, 1995); 3) the Indosinian Orogeny, South China (~60 million years from 270 to 210 Ma; e.g. Gilder et al., 1996). Furthermore, Duclaux et al. (2007) reported an integration study between triaxial numerical experiments and Archean field examples, and showed that strain rates in Archaean continental lithosphere might be in the range from ~5 × 10−15 to 10−16 s−1, allowing the preservation of L-S fabrics. However, there is a lack of any reliable data about true strain rates for the Yilgarn. The lower strain rates (0.5 × 10−15 and 1.0 × 10−15 s−1) considered in the present models reflect the lower bound for an active plate margin as stated above in the model boundary condition section. The modelling results show that the domination of strain localisation around the hot intrusion at the strain rates of ~10−15 s−1 only occurred at the very early stages of deformation (Figs. 2a and 4a). Although these sites of early strain localisation affect further strain localisation (i.e. acting as nucleation points for shear zones), the development of broad crustal scale shear zones in the models are not confined to the domains or boundaries of the intrusion (see Figs. 2c, d, 4c, d, 6c, d, g and h). These shear zones are well developed into the “cold” crustal domains (greenstone and granitic units) and they utilize lithological contacts (e.g. the flanks of the “cold” granite dome). This is due to the fact that the cooling of the intrusion to background temperatures occurred over a much shorter time interval when compared to the duration of the deformation event and development of shear zones. The cooling path of the hot intrusion for an extensional model with a strain rate of 1.0 × 10−15 s−1 (Fig. 8) shows that the intrusion was already cooled by more than 300 °C near the background geothermal state by 1% bulk extension. This means that for the main deformation period of the model, the intrusion had already cooled, with only some low remnant intrusive heat remaining (see Figs. 3 and 5).

Fig. 8. Plot of temperature versus time variations for a location within the magmatic intrusion, illustrating the cooling path of the intrusion. The model has extensional boundary conditions and a strain rate of 1.0 × 10−15 s−1 (see Figs. 2 and 3 for the corresponding strain and temperature distribution).

5.2. Strain localisation and implications on gold mineralisation An important observation from the current model results is that strong strain localisation near the intrusion and the development of shear zones in the crust occurred under both extensional and compressional tectonic conditions. There is always clear strain localisation around the shoulder areas of the intrusion and along the flanks of the “cold” granitic dome in the very early stages of deformation (Figs. 2a, b and 3a, b). Further structural development is dominated by shear zone development through the crust, which structurally utilizes or grows from earlier higher strain spots. For most strain rate cases, the final architecture is a graben basin for the extensional scenario (Figs. 2d and 6b–d) and a thickened “pop-up” wedge structure (Figs. 4d and 6f–h) for the compressional scenario, where these structures are bounded by normal-sense and reverse-sense shear zones/thrusts, respectively. Such graben basins are broadly consistent with the geometries of fault grabens in extensional basins from seismic interpretation (e.g. Gartrell et al., 2005) or experimental work (e.g. Dubois et al., 2002), and the “pop-up” structure also mimics the features of compressional positive flower-like “popup” structures due to reverse faulting and thrusting (e.g. McClay and Bonora, 2001; Viola et al., 2004) and associated high topography development (e.g. Huang et al., 2009). The extreme topography created at the top of these models reflects an over-estimated amplitude, which would be reduced if erosion had been considered in the model. The strain localisation features above provide a generic explanation for the observations of structural controls on the development of certain Au-mineralised sites in the Yilgarn Craton (e.g. Blewett et al., 2010a). Gold mineralisation in the Yilgarn occurred from the D2 contractional event (Tarmoola Gold deposit) and major gold mineralisation started from the D3 extensional event (e.g. the Sons of Gwalia deposit) through to the D4 contractional-transpressional events (e.g. New Holand deposit). During this series of extensional and contractional alternative deformation events, the emplacement of magmatic intrusions occurred approximately from D2/M2 (metamorphic) events onwards (e.g. Blewett et al., 2010a; Goscombe et al., 2007). This means that the emplacement of intrusions would be associated with either extensional or contractionaltranspressional deformation. Based on the current modelling results, such tectonic and magmatic conditions are favourable for strain localisation around the shoulder areas of hot intrusions and along the flanks of earlier granite domes. These strain localisation sites can be parts of broad crustal-scale shear zones which provide fluid conduits for mineralisation. In addition, the geometrical features here explain why gold deposits are structurally located at the flanks and apex of granite intrusion or dome (e.g. Blewett et al., 2010b; Henson et al., 2007; Henson and

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

Blewett, 2006). These areas are able to focus localisation over a protracted period and during changing deformation conditions. In the Yilgarn, fault structures were already developed before and during the contractional D2 deformation event (e.g. Blewett et al., 2010a). These faults became reactivated from the D3 events onwards, together with the development of other structures and gold mineralisation as described in Section 3.1. Therefore, the deformation–intrusion models with a pre-existing fault generically represent an important scenario for the Yilgarn. The results for the models with a pre-existing fault (Fig. 7) show that under both extensional and compressional conditions, early asymmetrical strain localisation occurred around the shoulder areas of the intrusion, with much stronger strain localisation at one shoulder adjacent to the reactivated fault. This is expressed as the formation of a damage zone linking the intrusion to the reactivated fault. Such structural patterns suggest that the flank of the intrusion adjacent to a pre-existing and reactivated fault would represent more favourable structural sites for gold mineralisation. This explains in part the structural control on gold mineralisation in the Yilgarn, such as the Wallaby deposit (Henson et al., 2007, 2010) and Paddington deposit (Davis et al., 2010), which are located near the flank of granite-cored domes and deep “tapping” faults (also see Blewett et al., 2010b for the description of gold mineralisation in the footwall of master faults and near domes). The models also show that such structures can develop under both extensional and convergent geodynamic conditions, in agreement with the geodynamic nature of major gold phases in the Yilgarn. Other pertinent structural features in the extensional model (Fig. 7a–c), include the large normal movement along the reactivated fault which generated fault blocks and a half graben. This simulates the formation of late basins in the Yilgarn, which occurred adjacent to the hanging-wall of reactivated faults as half graben or fault graben during the D3 extensional events (Blewett et al., 2010a; also see Goscombe et al., 2007). For compressional conditions (Fig. 7d–f), the large reverse or thrusting movement of the reactivated fault led to the development of multiple “pop-up” wedge structures or fault blocks. It is particularly interesting to note that movements along two secondary shear zones (Fig. 7f) switched to normal shearing and resulted in the development of a small down-thrown fault block or fault basin within the broad “pop-up” fault block under an overall contractional setting. This might reflect one possible mechanism for the formation of faulted basins under compressional conditions. We also observed in both extensional and compressional models, that new shear zones were developed along with the reactivation of the pre-existing fault and bulk deformation, in a style consistent with the concepts for antithetic faults under extension (e.g. F1 in Fig. c) and back-thrust for compression (e.g. F2 in Fig. 7f). Previous, field-based structural studies of shear zones and their relationship to gold mineralisation (e.g. Miller et al., 2010; Weinberg et al., 2004, 2005; Williams et al., 1989) have established a significant understanding of the controls of shear zones on gold mineralisation in the Yilgarn. Fluid inclusion study for shear zones and gold deposits in the Yilgarn (e.g. Baker et al., 2010) also generated insights on the pressure-temperature-compositional conditions for mineralising fluids around shear zones. In addition, there are also attempts to generically model fluid flow and gold precipitation patterns around a fault/shear zone by deformation, fluid flow and chemical reaction models (e.g. Zhang et al., 2007a; Zhao et al., 2007). The results of this study do not make direct contribution to the understanding of fluid flow and gold mineralisation in shear zones. The contributions of the current modelling work are mainly reflected in two aspects: (1) the Yilgarn application models demonstrated that this new continuum-damage-mechanics algorithm and approach can realistically simulate the localisation of very high strain and the formation of shear zones; and (2) the strain localisation and shear zone patterns of the models generically explain why gold deposits in the Yilgarn occurred (observed) at certain structural locations as described above in a generic manner.

145

6. Conclusions This study has simulated rock deformation around a magmatic intrusion for a number of strain rate and structural scenarios relevant to the Yilgarn, employing a new continuum damage mechanics algorithm that considers the temperature and time dependent elastic-viscoplastic constitutive behaviour of crustal materials. The model results show that strain rates have the most significant effects on structural patterns for both extensional and compressional cases. Low strain rates (10−15 s −1) promote the formation of narrow high-strain shear zones and strong strain localisation along the flanks or shoulder areas of the intrusion and cold granite dome. Under higher strain rate conditions, diffuse shear zones (wider and lower strains) developed as the result of strain hardening. The cooling rate of the intrusion is more rapid than bulk deformation. At a strain rate of 1.0 × 10 −15 s−1 (a rate applicable to many orogenic events), for example, the intrusion was cooled by more than 300° within 1% bulk extension, suggesting that for the main deformation period, the intrusion would become essentially “cooled” near the background geothermal state with only some low remnant intrusive temperature. Strong strain localisation near the intrusion and the development of shear zones in the crust occurred under both extensional and compressional conditions. There is always a clear strain localisation around the shoulder areas of the intrusion and along the flanks of the “cold” granitic dome in the very early stages of deformation. These strain localisation sites affected shear zone development through the crust in later stages. In the cases where a pre-existing fault was included, early strain localisation near the intrusion became asymmetrical and there is much stronger strain localisation at one shoulder adjacent to the reactivated fault for both extensional and compressional settings. This is shown as the formation of a damage zone linking one flank of the intrusion to the reactivated fault, suggesting more favourable structural conditions for gold mineralisation there. At higher bulk deformation stages, the models produced the structural patterns of graben basin (extensional scenario) and “pop-up” wedge structure (compressional scenario), bounded by normal and reverse shear zones, respectively. Large normal movement along the reactivated fault in the extensional model generated the geometries of tilted fault blocks and a half graben, which represents a possible mechanism for the formation of late basins in the Yilgarn. The results of the compressional model show that switch from reverse to normal movement and development of local faulted basins are possible within a “pop-up” structure under broad compressional conditions. The model results provide insights on the development of strain localisation and shear zones at certain structurally-favourable locations in the Yilgarn, which are associated with gold mineralisation. This may explain why a number of gold deposits in the Yilgarn (e.g. the Wallaby deposit and Paddington deposit) were formed near the flank of granite-cored domes and deep “tapping” faults. The models also show that these structures can occur under both extensional and convergent tectonic conditions, consistent with the geodynamic features of major gold mineralisation phases reported for the Yilgarn Craton. The consistency between the model results and certain geological observations shows that this new continuum-damagemechanics algorithm and modelling approach are capable of realistically simulating high strain localisation and shear zone development. Acknowledgement We acknowledge the funding support from the CSIRO-led Minerals Down Under National Research Flagship on this study. References Ashby, M.F., Gandhi, C., Taplin, D.M.R., 1979. Overview No. 3 Fracture-mechanism maps and their construction for f.c.c. metals and alloys. Acta Metallurgica 27, 699–729.

146

Y. Zhang et al. / Tectonophysics 526–529 (2012) 133–146

Baker, T., Bertelli, M., Blenkinsop, T., Cleverley, J.S., McLellan, J., Nugus, M., Gillen, D., 2010. P–T–X conditions of fluids in the Sunrise Dam Gold Deposit, Western Australia, and implications for the interplay between deformation and fluids. Economic Geology 105, 873–894. Blewett, R.S., Czarnota, K., Henson, P.A., 2010a. Structural-event framework for the eastern Yilgarn Craton, Western Australia, and its implications for orogenic gold. Precambrian Research 183, 203–329. Blewett, R.S., Henson, P.A., Roy, I., Champion, D.C., Cassidy, K.F., 2010b. Scale-integrated architecture of a world-class gold mineral system: the Archaean eastern Yilgarn Craton, Western Australia. Precambrian Research 183, 230–250. Cassidy, K.F., Hagemann, S.G., 2001. ‘World-class’ Archean orogenic gold deposits, eastern Yilgarn Craton: diversity in timing, structural controls and mineralization styles. AGSO-Geoscience Australia. Record 2001 (37), 382–384. Cassidy, K.F., Groves, D.I., McNaughton, N.J., 1998. Late-Archean granitoid-hosted lodegold deposits, Yilgarn Craton, Western Australia: Deposit characteristics, crustal architecture and implications for ore genesis. Ore Geology Reviews 13, 65–102. Czarnota, K., Blewett, R.S., Goscombe, B., 2010a. Predictive mineral discovery in the eastern Yilgarn Craton, Western Australia: an example of district scale targeting of an orogenic gold mineral system. Precambrian Research 183, 356–377. Czarnota, K., Champion, D.C., Goscombe, B., Blewett, R.S., Cassidy, K.F., Henson, P.A., Groenewald, P.B., 2010b. Geodynamics of the eastern Yilgarn Craton. Precambrian Research 183, 175–202. Davis, B.K., Blewett, R.S., Squire, R., Champion, D.C., Henson, P., 2010. Granite-cored domes and gold mineralisation: Architectural and geodynamic controls around the Archaean Scotia-Kanowna Dome, Kalgoorlie Terrane, Western Australia. Precambrian Research 183, 316–337. Dubois, A., Odonne, F., Massonnat, G., Lebourg, T., Fabre, R., 2002. Analogue modelling of fault reactivation: tectonic inversion and oblique remobilisation of grabens. Journal of Structural Geology 24, 1741–1752. Duclaux, G., Rey, P., Guillot, S., Menot, R.-P., 2007. Orogen-parallel flow during continental convergence: Numerical experiments and Archean field examples. Geology 35, 715–718. Dunster, J.N., Kruse, P.D., Duffett, M.L., Ambrose, G.J., 2007. Geology and resource potential of the southern Georgina Basin. Northern Territory Geological Survey, Report DIP 007. 232 pp. Fusseis, F., Regenauer-Lieb, K., Liu, J., Hough, R., deCarlo, F., 2009. Creep Cavitation can establish a granular fluid pump through the middle crust. Nature 459, 974–977. Gartrell, A., Bailey, W., Brincat, M., 2005. Strain localisation and trap geometry as key controls on hydrocarbon preservation in the Laminaria High area. APPEA Journal 45, 477–492. Garven, G., Bull, S.W., Large, R.R., 2001. Hydrothermal fluid flow models of stratiform ore genesis in the McArthur Basin, Northern Territory, Australia. Geofluids 1, 289–311. Gilder, S.A., Gill, J., Coe, R.S., Zhao, X.X., Liu, Z.W., Wang, G.X., 1996. Isotopic and paleomagnetic constraints on the Mesozoic tectonic evolution of South China. Journal of Geophysical Research 107 (B7), 16137–16154. Goldfarb, R.J., Groves, D.I., Gardoll, S., 2001. Orogenic gold and geologic time: a global synthesis. Ore Geology Reviews 18, 1–75. Goscombe, B., Blewett, R., Czarnota, K., Maas, R., Groenewald, B., 2007. Broad thermobarometeric evolution of the Eastern Goldfields Superterrane. In: Bierlein, F.P., Knox-Robinson, C.M. (Eds.), Proceedings of Geoconferences (WA) Inc. Kalgoorlie'07 conference: Geoscience Australia Record 207/14, pp. 33–38. Groves, D.I., Goldfarb, R.J., Know-Robinson, C.M., Ojala, J., Gardoll, S., Yun, G.Y., Holyland, P., 2000. Late-kinematic timing of orogenic gold deposits and significance for computer-based exploration techniques with emphasis on the Yilgarn Block, Western Australia. Ore Geology Reviews 17, 1–38. Henson, P.A., Blewett, R.S., 2006. Question 5 – metal transport and depositional processes. In: Blewett, R.S., Hitchman, A.P. (Eds.), Final report, 3D geological models of the eastern Yilgarn Craton, Project Y2: Geoscience Australia Record 2006/5, pp. 242–252http://www.ga.gov.au/images_cache/GA13801.pdf. Henson, P.A., Blewett, R.S., Champion, D.C., Goleby, B.R., Czarnota, K., 2007. How does the 3D architecture of the Yilgarn control hydrothermal fluid focusing? In: Bierlein, F.P., Knox-Robinson, C.M. (Eds.), Proceedings of Geoconferences (WA) Inc Kalgoorlie'07 conference: Geoscience Australia Record 207/14, pp. 57–61. Henson, P.A., Blewett, R.S., Roy, I.G., Miller, J.M., Czarnota, K., 2010. 4D architecture and tectonic evolution of the Laverton region, eastern Yilgarn Craton, Western Australia. Precambrian Research 183, 338–355. Hobbs, B.E.H., Zhang, Y., Ord, A., Zhao, C., 2000. Application of coupled deformation, fluid flow, thermal and chemical modelling to predictive mineral exploration. Journal of Geochemical Exploration 69–70, 505–509. Holyland, P.W., Ojala, V.j., 1997. Computer aided structural targeting in mineral exploration: two to three-dimensional stress mapping. Australian Journal of Earth Sciences 44, 421–432. Huang, M.H., Hu, J.C., Ching, K.e., Rau, R.J., Hsieh, C.S., Pathier, E., Fruneau, B., Deffontaines, B., 2009. Active deformation of Tainan tableland of southwestern Taiwan based on geodetic measurements and SAR interferometry. Tectonophysics 466, 322–334. Ji, S., Zhao, P., Xia, B., 2003. Flow laws of multiple materials and rocks from end-member flow laws. Tectonophysics 370, 129–145. Karrech, A., Regenauer-Lieb, K., Poulet, T., 2011a. Continuum damage mechanics for the lithosphere. Journal of Geophysical Research B: Solid Earth 116 (4), 1–14 art. no. B04205. Karrech, A., Regenauer-Lieb, K., Poulet, T., 2011b. Frame indifferent elastoplasticity of frictional materials at finite strain. International Journal of Solids and Structures 48 (3–4), 397–407. Karrech, A., Regenauer-Lieb, K., Poulet, T., 2011c. A damaged visco-plasticity model for pressure and temperature sensitive geomaterials. International Journal of Engineering Science 49, 1141–1150.

Kreemer, C., Holt, W.E., Haines, A.J., 2003. An integrated global model of present-day plate motions and plate boundary deformation. Geophysical Journal International 154, 8–34. Lee, T.Y., Lawver, L.A., 1995. Cenozoic plate reconstruction of Southeast Asia. Tectonophysics 251, 85–138. McClay, K., Bonora, M., 2001. Analog models of restraining stepovers in strike-slip fault systems. AAPG Bulletin 85, 233–260. Miller, J., Blewett, R.S., Tunjic, J., Connors, K., 2010. The role of early formed structures on the development of the world class St Ives Goldfield, Yilgarn, WA. Precambrian Research 183, 292–315. Ord, A., Hobbs, B.E., Zhang, Y., Broadbent, G.C., Brown, M., Willetts, G., Sorjonen-Ward, P., Walshe, J.L., Zhao, C., 2002. Geodynamic modeling of the Century deposit, Mt Isa Province, Queensland. Australian Journal of Earth Sciences 49, 1011–1039. Potma, W., Roberts, P.A., Schaubs, P.M., Sheldon, H.A., Zhang, Y., Hobbs, B.E., Ord, A., 2008. Predictive targeting in Australian orogenic-gold systems at the deposit to district scale using numerical modelling. Aust J Earth Sciences 55, 101–122. Ranalli, G., 1997. Rheology and deep tectonics. Annali di Geofisica XL (3), 671–680. Ranalli, G., 2000. Rheology of the crust and its role in tectonic reactivation. Journal of Geodynamics 30, 3–15. Regenauer-Lieb, K., 1999. Dilatant plasticity applied to Alpine Collision: ductile void growth in the intraplate area beneath the Eifel Volcanic Field. Journal of Geodynamics 27, 1–21. Regenauer-Lieb, K., Weinberg, R.F., Rosenbaum, G., 2006. The effect of energy feedbacks on continental strength. Nature 442, 67–70. Regenauer-Lieb, K., Karrech, A., Chua, H.T., Horowitz, F.G., Yuen, D., 2010. Time-dependent, irreversible entropy production and geodynamics, Philosophical Transactions of the Royal Society A: Mathematical. Physical and Engineering Sciences 368 (1910), 285–300. Schaubs, P.M., Rawling, T.J., Dugdale, L.J., Wilson, C.J.L., 2006. Factors controlling the location of gold mineralisation around basalt domes in the Stawell corridor: insights from coupled 3D deformation-fluid-flow numerical models. Australian Journal of Earth Sciences 53, 841–862. Sorjonen-Ward, P., Zhang, Y., Zhao, C., 2002. Numerical modeling of orogenic processes and gold mineralization in the southeastern part of the Yilgarn craton, Western Australia. Australian Journal of Earth Sciences 49, 1011–1039. Turcotte, D.L., Schubert, G., 1982. Geodynamics: applications of continuum physics to geological problems. Wiley, New York. 450 pp. Viola, G., Odonne, F., Mancktelow, N.S., 2004. Analogue modelling of reverse fault reactivaton in strike-slip and transpressive regimes: application to the Giudicarie fault system, Italian Eastern Alps. Journal of Structural Geology 36, 401–418. Wang, Y.J., Fan, W.M., Cawood, P.A., Ji, S.C., Peng, T.P., Chen, X.Y., 2007. Indosinian high strain deformation for the Yunkaidashan tectonic belt, South China: kinematics and 40Ar/39Ar geochronological constraints. Tectonics 26 (TC6008). doi:10.1029/ 2007, TC002099. Weinberg, R., Regenauer-Lieb, K., 2010. Ductile fracture and magma migration from source. Geology 38, 363–366. Weinberg, R.F., van der Borgh, P., 2008. Extension and gold mineralization in the Archean Kalgoorlie Terrane, Yilgarn Craton. Precambrian Research 161, 77–88. Weinberg, R.F., Hodkiewicz, P.F., Groves, D.I., 2004. What controls gold distribution in Archean terranes? Geology 32, 545–548. Weinberg, R.F., van der Borgh, P., Bateman, R.J., Groves, D.I., 2005. Kinematic history of the Boulder-Lefroy shear zone system and controls on associated gold mineralization, Yilgarn Craton, Western Australia. Economic Geology 100, 1407–1426. Wibberley, C.A.J., Yielding, G., Di Toro, G., 2008. Recent advances in the understanding of fault zone internal structure: a review. In: Wibberley, C.A.J., Kurz, W., Imber, J., Holdsworth, R.E., Collettini, C. (Eds.), The internal structure of fault zones: implication for mechanical and fluid-flow properties. : Special Publications, 299. Geological Society, London, pp. 5–33. Williams, P.R., Nisbet, B.W., Etheridge, M.A., 1989. Shear zone, gold mineralization and structural history in the Leonora district, Eastern Goldfields Province, Western Australia. Australian Journal of Earth Sciences 36, 383–403. Zhang, Y., Hobbs, B.E., Ord, A., Barnicoat, A., Zhao, C., Walshe, J.L., Lin, Ge, 2003. The influence of faulting on host-rock permeability, fluid flow and mineral precipitation: a conceptual 2D numerical model. Journal of Geochemical Exploration 78–79, 279–284. Zhang, Y., Sorjonen-Ward, P., Ord, A., Southgate, P., 2006. Fluid flow during deformation associated with structural closure of the Isa Superbasin after 1575 Ma in the central and northern Lawn Hill Platform, Northern Australia. Economic Geology 101, 1293–1312. Zhang, Y., Cleverley, J.S., Sheldon, H.A., Ord, A., Blewett, R.S., Barnicoat, A.C., 2007a. Modeling of rock deformation, fluid flow, reactive transport and Au mineralization: examples relating to the Yilgarn craton, Australia. In: Bullen, T.D., Wang, Y. (Eds.), Water-Rock Interaction (Proceedings of the 12th International Symposium on Water-rock interaction, Kunming, China, 31 July-5 August 2007). Taylor & Francis Group, London, pp. 175–179. ISBN 978-00415-45136-9. Zhang, Y., Lin, G., Roberts, P.A., Ord, A., 2007b. Numerical modelling of deformation and fluid flow in the Shui-Kou-Shan mineralisation district, Hunan Province, South China. Ore Geology Review 31, 261–278. Zhao, C., Hobbs, B.E., Mühlhaus, H.B., 1998. Finite element modelling of temperature gradient driven rock alteration and mineralization in porous rock masses. Computer Methods in Applied Mechanics and Engineering 165, 175–187. Zhao, C., Hobbs, B.E., Ord, A., Hornby, P., Peng, S., Liu, L., 2007. Mineral precipitation associated with vertical fault zones: the interaction of solute advecton, diffusion and chemical kinetics. Geofluids 7, 3–18.