- Email: [email protected]
Three-dimensional numerical modeling of salinity variations in driving basin-scale ore-forming fluid flow: Example from Mount Isa Basin, northern Australia
Journal of Geochemical Exploration 106 (2010) 236–243
Contents lists available at ScienceDirect
Journal of Geochemical Exploration j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j g e o ex p
Three-dimensional numerical modeling of salinity variations in driving basin-scale ore-forming fluid flow: Example from Mount Isa Basin, northern Australia Jianwen Yang a,b,⁎, Zuohai Feng a, Xianrong Luo a, Yuanrong Chen a a b
Sabbatical leave at Faculty of Earth Sciences, Guilin University of Technology, Guilin, Guangxi, 541004, PR China Department of Earth and Environmental Sciences, University of Windsor, Windsor, Ontario, Canada N9B 3P4
a r t i c l e
i n f o
Article history: Received 27 April 2009 Accepted 16 December 2009 Available online 28 December 2009 Keywords: Hydrothermal fluid flow Finite element modeling SEDEX deposits Mount Isa basin
a b s t r a c t This paper presents a fully 3-D numerical investigation into the effect of salinity on ore-forming hydrothermal fluid flow, and develops a highly conceptualized hydrological model to simulate the fluid plumbing system early in the history of the Mount Isa basin, northern Australia when lead–zinc deposits were formed therein. Our numerical modeling results indicate that the active synsedimentary faults and clastic aquifer form a favourable hydro-framework for regional-scale fluid flow, and variations in salinity have important implications for fluid migration and heat transport. When salinity is constant throughout the basin, hydrothermal fluid flow mainly circulates within the more permeable faults as 2-D convection cells, unless the contrast in permeability between the faults and aquifer is less than one order of magnitude. Enhanced salinities on the basin floor due to evaporation of seawater facilitate the development of full 3-D thermohaline flow systems. The 3-D convection rolls established by the evaporitic conditions tend to be stretched if only longitudinal faults exist, but become ‘mushroom-shaped’ when a cross fault is added to intersect the longitudinal ones. The intersection of the two sets of the faults enables seawater to channel downwards to a greater depth and circulate through larger volumes of the source rocks to leach more metal content, and eventually the heated basinal brines ascend to the basin floor via a more ‘localized’ discharge zone at higher venting fluid velocity and temperature, potentially forming giant ore deposits. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Salinity (evaporitic) conditions of sedimentation have long been recognized as an important factor in the development of basinal brines that form stratiform Zn–Pb–Ag ore deposits, mainly because the solubility of the economically important metals increases with increasing salinity (e.g., Hanor, 1996). Also, spatial and temporal changes in salinity can lead to a buoyancy force due to density variation of basinal brine, which is critical in determining regional fluid flow patterns. Buoyancy-driven fluid flow has been extensively studied analytically and numerically in large-scale geological systems (e.g., Zhao et al., 1997, 1999). It is recently gaining momentum as the most likely hydrological scenario for the formation of sedimentary-exhalative (SEDEX) ore deposits. Previous numerical studies usually considered thermallyinduced fluid flow only, with no account taken of the effects of salinity (e.g., Barrier et al., 1999; Garven et al., 2001; Yang et al., 2006). A few studies have indeed considered both the heat- and salinity-induced buoyancy effects on fluid flow in porous media that are heated and ⁎ Corresponding author. Department of Earth and Environmental Sciences, University of Windsor, Windsor, Ontario, Canada N9B 3P4. Tel.: +1 519 253 3000x2181; fax: +1 519 973 7081. E-mail address: [email protected] (J. Yang). 0375-6742/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.gexplo.2009.12.004
salted from below due to the co-presence of natural geothermal gradients and salt domes (e.g., Sarkar et al., 1995; Schoofs and Spera, 2003), but they are not related to SEDEX ore genesis. We have recently developed a hydrogeological model that fully couples transient fluid flow, heat and solute transport associated with the formation of the HYC SEDEX deposit in the McArthur Basin, northern Australia (Yang et al., 2004; Yang, 2006a). However, these numerical experiments are 2-D in nature over the selected geological sections, and therefore cannot fully represent the complexities of 3-D hydrothermal system in reality. Recent fast developments in computer hardware and software have enabled researchers to start with simulating hydrothermal ore-forming fluid flow in realistic 3-D environments. For instance, Zhao et al. (2003) employed the finite element model to simulate 3-D steady-state convective pore-fluid flow and the relevant mineralization in fluidsaturated rocks. Zhao et al. (2006) derived analytical and numerical solutions for the double-diffusion driven convective instability in 3-D permeable faults, once again under steady-state conditions. Yang (2006b) presented a 3-D hydrological model that fully couples transient fluid flow and heat transport in the McArthur basin, northern Australia, but without considering the effect of salinity distribution on regional fluid flow. Zhao et al. (2008a) employed numerical modeling to address the controlling formation mechanism related to equal-distant distributed gold deposits in a 3-D fault. Zhao et al. (2008b) considered the morphological evolution of 3-D chemical dissolution fronts occurring in
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
fluid-saturated porous media. More recently Feltrin et al. (2009) simulated 3-D fluid flow driven by tectonic deformation associated with the formation of the giant lead–zinc–silver Century deposit, but no account was taken of the effect of temperature and salinity variation on fluid flow in the hydrothermal system. This paper presents a 3-D numerical investigation that fully couples transient fluid flow and heat transfer with solute transport, and it concerns in particular with buoyancy-driven fluid flow pertinent to the genesis of shale-hosted lead–zinc ores in the Mount Isa basin, northern Australia. In this study, we do not examine any specific deposits; rather we develop a highly generalized 3-D conceptual model on the basis of the general geological structure of the targeted basin. In addition in order to simplify the physics involved, we only consider a single solute transport and ignore chemical reactions between different phases. Ore-forming fluid flow can be driven by tectonic deformation as well (e.g., Zhang et al., 2009), however we do not intend to address this issue as our attention is focused on highlighting the role of salinity variations in driving basinal brines. We conduct a series of numerical experiments with different conditions in order to address the effect of salinity variations in driving hydrothermal ore-forming fluid migration in the study region. 2. Governing equations and finite element modeling Mathematically, saline hydrothermal fluid flow in sedimentary basins is governed by the fluid continuity equation, thermal energy conservation equation, and solute mass conservation equation, coupled through Darcy's law. The fluid continuity equation can be expressed in terms of an ‘equivalent freshwater’ hydraulic head for the variable density fluid flow system (Yang, 2006a). For a 3-D system in the x–y–z coordinate, we have: ∂ ∂h ∂ ∂h ∂ ∂h ∂h K + K + K + Kρr = Ss ; ∂x ∂x ∂y ∂y ∂z ∂z ∂t ∂ ∂T ∂ ∂T ∂ ∂T λm + λm + λm ∂x ∂x ∂y ∂y ∂z ∂z −
∂ ∂ ∂ ∂T ðc ρ q TÞ− ðcw ρw qx TÞ− ðcw ρw qz TÞ = cm ρm ; ∂x w w x ∂y ∂z ∂t
ð1Þ
ð2Þ
and ∂ ∂C ∂ ∂C ∂ ∂C θDx + θDy + θDz ∂x ∂x ∂y ∂y ∂z ∂z −
∂ ∂ ∂ ∂ ðq CÞ− ðqy CÞ− ðqz CÞ = ðθRCÞ; ∂x x ∂y ∂z ∂t
ð3Þ
where qx, qy and qz are the Darcy flux components in x, y, and z directions, h is the ‘equivalent freshwater’ head, K is the hydraulic − θ) conductivity, Ss is the specific storage, λm = λθwλ(1 (λw and λs s denote the thermal conductivity of the fluid and solid phase), θ is the porosity, cmρm = cwρwθ + csρs(1 − θ) (cw and cs are the specific heat capacity of the fluid and solid phase, and ρw and ρs are the density of the fluid and solid phase), ρr is the relative fluid density defined as ρr = ρw / ρ0 − 1 (ρ0 is the reference fresh water density), t is the time, C is the solute concentration, T is the temperature, R is the retardation factor, and Dx, Dy and Dz are the hydrodynamic dispersion coefficients. The governing Eqs. (1)–(3) are not alone sufficient. The dependence of fluid density and viscosity on temperature and solute concentration must be defined. In this paper it is calculated using the NIST/ASME Steam Properties code (Klein and Harvey, 1996). These equations form a time-dependent, nonlinear and coupled system, which leads their solutions to become nontrivial even for a very simple 1-D geological system. We have recently developed a finite element computer model, following original works by Molson et al. (1992), to numerically solve these governing equations. Listed
237
here are the major modeling steps: 1) update the fluid density and viscosity using the latest solute concentration and temperature from the previous time step; 2) determine the ‘equivalent freshwater’ head based on the updated fluid properties by solving the fluid continuity equation; 3) determine the Darcy flux based on Darcy's law using the updated fluid properties and hydraulic head; 4) calculate the temperature and solute concentration by respectively solving the thermal energy and solute mass conservation equations; and 5) check whether the preset convergence criteria are satisfactory. If not, steps 1) through 4) are repeated; if yes, move on to the next time step. We adopt a non-orthogonal quadrilateral mesh to create finite elements since it is better suited for the complex geometry of different stratigraphic units and faults encountered in sedimentary basins. Each element is assigned permeability, porosity, thermal conductivity and other physical parameters governing fluid flow, heat transport and solute transport based on the rock properties. Further details of the computational method, both the principles and the schemes for their implementation into finite element software, can be found from our previous publications (Yang et al., 2004; Yang, 2006a). 3. Three-dimensional conceptualized model As illustrated in Fig. 1a, our 3-D conceptual model has dimensions of 60, 30, and 20 km respectively in the x, y, and z directions, and is discretised by a 3-D non-orthogonal quadrilateral mesh consisting of 21,736 finite elements in total. Fig. 1b illustrates the paleohydrostratigraphy over the central x–z cross-section, which is constrained by some of the common features of a sedimentary basin's rift-and-sag phase, and in particular by the reconstructions of the Mount Isa basin, northern Australia (O'dea et al., 1997; Betts et al., 2003). Fig. 1 represents the highly conceptualized and simplified subsurface stratigraphy and structure, which controlled the hydrological system when the Mount Isa SEDEX deposits were formed early in the history of the Mount Isa basin. As shown in Fig. 1b, the conceptual model involves a volcanic basement sequence of low permeability (Unit 1), a sandstone aquifer of high permeability (Unit 2), a rift cover sequence of intermediate permeability (Unit 3), and an upper cover sequence of shales and siltstones (Unit 4) that hosts mineral deposits formed during or soon after sediment deposition. Two more permeable faults (Fault 1 and Fault 2) are also included to penetrate from the upper sequence into the basement. The faults are assumed to be 1 km wide and steeply-dipping, which is constrained by surface exposure and seismic profiling (Bierlein and Betts, 2004). A similar general geological structure of the Mount Isa–McArthur basin region has been also used in previous fluid flow models (e.g., Yang et al., 2004; McLellan et al., 2006; Oliver et al., 2006; Yang, 2006a). Like in the previous studies of McLellan et al. (2006) and Oliver et al. (2006), Fault 1 on the left is equivalent to the Mount Isa fault system, and Fault 2 on the right simulates a fault zone to farther north, such as the Termite Range fault at the Century deposit. The two steeplydipping faults cut the sandstone aquifer (Unit 2), forming a favourable hydrological framework for regional-scale fluid flow in both longitudinal and transverse directions. The permeabilities and thermal conductivities assigned to the stratigraphic elements and faults are given in Table 1, following the previous numerical investigations in the Mount Isa–McArthur basin region listed above. We assume that the vertical permeability of the host rocks is two orders of magnitude less than the horizontal permeability due to the bedded and stratified nature of these sedimentary rocks. In addition, we assume cw = 4174 J/kg °C (specific heat of the fluid phase), cs = 800 J/kg °C (specific heat of the solid phase), λw = 0.5 W/m °C (thermal conductivity of the fluid phase), ρs = 2630 kg/m3 (density of the solid phase), ρ0 = 1000 kg/m3 (reference fresh water density), and θ = 10% (porosity). The upper boundary (i.e., the basin floor) is permeable to fluid flow and its temperature is fixed at 20 °C. The lower boundary temperature is maintained at 450 °C, justified by present-day heat
238
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
Fig. 1. 3-D conceptualized hydrological model based on the general geological structure of the Mount Isa basin, northern Australia: (a) finite element mesh, and (b) central crosssection. Refer to Table 1 for the hydrological and thermal properties assigned to the faults and rock units.
flow measurements in northern Australia (S. McLaren and M. Sandiford, pers. comm.). The bottom of the model is assumed impermeable since it lies within the deep volcanic basement. The side boundaries are assumed to be adiabatic to heat transfer and impermeable to fluid flow. The initial fluid velocity is set to zero over the whole solution domain, and the initial temperature is assumed to vary linearly with depth. Various salinity conditions are considered in our numerical experiments as follows.
4. Numerical simulations and results The first scenario corresponds to a submarine condition, as geological evidence (e.g., Yang et al., 2006) indicates that open marine conditions persisted during deposition of the upper Isa Superbasin. We assume that basinal brine in this case is simply pure seawater, and therefore assign a constant salinity of typical seawater Table 1 Major physical parameters of the faults and host rocks. Hydrological unit Horizontal permeability Vertical permeability Thermal conductivity and formation (m2) (m2) (J/m/s/°C) Fault Unit 1 Unit 2 Unit 3 Unit 4
4 × 10− 14 2 × 10− 16 4 × 10− 14 2 × 10− 15 1 × 10− 15
4 × 10− 14 2 × 10− 18 4 × 10− 16 2 × 10− 17 1 × 10− 17
2.0 3.0 3.0 2.5 2.0
3.57% throughout the basin. Clearly under this condition, fluid flow is driven by thermally-induced buoyancy force only. Fig. 2 shows the numerical results corresponding to the reference conditions specified above and given in Fig. 1 and Table 1. Fig. 2a and b illustrates the temperature evolution from 0.5 to 0.8 Ma, and Fig. 2c shows the fluid velocity vectors at the time of 0.8 Ma. It can be seen that due to the large contrast of hydraulic conductivities between the faults and host rocks, significant fluid circulation is confined mainly within the more permeable faults with a maximum rate of 2.3 m/year rather than in the host media (Fig. 2c). Both upwelling and downwelling flow develops within these two faults at different longitudinal (y-axis) distances, which gives rise to significant temperature variation along the fault strike direction (Fig. 2b). By comparison, fluid flow in the host rocks including the aquifer is less significant, evidenced by the nearly flat isothermal lines (Fig. 2a and b). This is very similar to our previous 3-D modeling results of the McArthur basin under pure water conditions (Yang, 2006b). Since hydrothermal fluid flow more easily focuses and circulates within the fault zones with less fluid transport in the host rocks, the metal contents leached from the host rocks by the circulating basinal brines are likely insufficient to form giant SEDEX deposits. However, this situation may change if the sandstone aquifer (Unit 2) is more permeable. In the second scenario, therefore, we assume that the vertical permeability of the aquifer is increased by one order of magnitude (but still one order of magnitude lower than the permeability of the faults). Fig. 3 shows the temperature distribution and fluid velocity vectors at the time of 0.8 Ma, from which it can be
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
239
Fig. 2. Scenario 1 corresponds to a constant salinity of typical seawater of 3.57% assigned throughout the basin: (a) temperature contours at a time of 0.5 Ma; (b) temperature contours at a time of 0.8 Ma; and (c) fluid velocity vectors over the faults at a time of 0.8 Ma.
seen that two 3-D convection rolls now develop. Cold seawater recharges the basin via the central segment of Fault 1 (Fig. 3b and c), and it then not only circulates longitudinally along the fault like in the first scenario but also travels laterally through the aquifer toward the right (Fig. 3a and c). In the meantime, the basinal fluid is heated up from below and eventually ascends via the central segment of Fault 2 and discharges to the basin floor to form a SEDEX-type deposit (Fig. 3b and c). Fluid discharge temperatures at the surface of the model are in the range of 135 °C to 190 °C with fluid velocities of 1.5 to 2.5 m/year
for the period of 0.8 Ma. Note that as a 3-D convection roll, cold seawater also recharges the basin via the off-central parts of Fault 2 (Fig. 3a and b). It circulates downwards via the fault and the surrounding volume of the host rocks, heated up at depth, and finally flows upwards and discharges to the basin floor via the central part of Fault 2 (Fig. 3b and c). Our further numerical experiments (results not shown here) indicate that when the permeability contrast between the faults and aquifer is less than one order of magnitude, 3-D convective fluid flow system can always develop to circulate over a
240
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
Fig. 3. Scenario 2 corresponds to a constant salinity of typical seawater of 3.57% assigned throughout the basin, and the vertical permeability of the aquifer Unit 2 is increased by one order of magnitude from the reference value in Table 1: (a) temperature contours at a time of 0.8 Ma; (b) iso-surface of temperature distribution at a time of 0.8 Ma; and (c) fluid velocity vectors over the central cross-section at the same time. Note that D and R denote the major discharge and recharge zones.
large volume of the host rocks for leaching metals from the source rocks, which is in favour of the ore genesis of SEDEX-type deposits. However, when the permeability contrast exceeds one order of magnitude, fluid flow tends to be confined within the more permeable fault zone as 2-D convection cells, and this is not favourable for the ore formation. The third scenario corresponds to semi-evaporitic conditions resulting from the mixing of minimally modified seawater in the
sub-basins with saline brines generated by evaporation at the basin margins. We double the salinity of typical seawater 3.57% for the upper boundary (basin floor), and assign the seawater salinity for the lower boundary. Initial salinity is also equal to 3.57% throughout the basin. Fig. 4a and b shows the snapshots of iso-surface of the salinities 6% and 7% at times of 0.05 Ma and 0.25 Ma, corresponding to the reference conditions given in Fig. 1 and Table 1. It can be seen that this type of saline conditions facilitates buoyancy-driven free convection,
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
241
Fig. 4. Scenario 3 corresponds to the semi-evaporitic condition: (a) iso-surface of salinity distribution at a time of 0.05 Ma; (b) iso-surface of salinity distribution at a time of 0.25 Ma; and (c) iso-surface of temperature distribution at a time of 0.25 Ma. Note that D and R denote the major discharge and recharge zones.
compared with the pure seawater case discussed above. Due to the continuous supply of brine from the seafloor and hence the enhanced convection strength, saline fluid circulates not only within the faults but also spreads out via the aquifer and host rocks, and as a result the salinity in the aquifer increases progressively with time. Consequently several 3-D convection rolls develop, as identified from the iso-surface of temperatures 80 °C and 240 °C at the time of 0.25 Ma shown in Fig. 4c, which allows the hot basinal brines to circulate through sufficient volumes of the host rocks to leach enough metal content.
Also, as illustrated in Fig. 4c, several elongated discharge zones develop on the top of the faults, over which SEDEX-type deposits may be formed. The venting fluid temperatures on the basin floor range from 180 °C to 210 °C with fluid velocities of 47 to 75 m/year over the period of 0.25 Ma, both greater than those shown in the above scenarios under the pure seawater condition. In the forth scenario, we include a vertical cross fault to reflect the overall effect of several minor faults that sinuously and intermittently cut the longitudinal major faults in the targeted region (e.g., refer to
242
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
Fig. 1 of Yang et al., 2006). Other conditions remain the same as in the third scenario. The cross fault, located in the central x–z plane at y = 15 km, orthogonally intersects the two major faults. To simplify the complexity, this cross fault is assumed to have the same width and hydraulic conductivity as those of the primary main faults. The intersection of these two sets of faults may result in localized zones of enhanced permeability; therefore, two subvertical columns, centered by the left-hand and right-hand intersections and having a crosssectional area of 3 km by 3 km, are included and assigned the same permeability of the faults. Fig. 5 illustrates the temperature distribution at the time of 0.25 Ma, from which it can be seen that fluids are mainly channelled downwards via the left-hand localized column and upwards through the right-hand conduit. This leads to the development of two major 3-D ‘mushroom-shaped’ convection rolls that are centered by these two columns of enhanced permeability, and two ‘minor’ 3-D convection rolls that are located close to the ends of the cross faults (Fig. 5b). Significant hydrothermal convection now takes place not only along the longitudinal direction as shown in Fig. 2 of the first case study, but also along the transverse direction through both the cross fault and the host rocks, as evidenced by several upwelling and downwelling thermal plumes distributed alternately over the cross fault (Fig. 5a). The F1-associated left intersection zone behaves as the major recharge pathway, whereas the F2-related right localized zone serves as the major discharge conduit. Cold seawater moves downwards along the left subvertical column, then travels laterally to the right through the cross fault and the aquifer (the host rock Unit 2)
and in the meantime it is heated from below, and eventually the hot basinal brines ascend via the right subvertical column and discharge to the basin floor to form a SEDEX-type deposit. The venting fluid temperatures on the basin floor range from 190 °C to 235 °C with fluid velocities of 66 to 138 m/year over the period of 0.25 Ma. Numerical modeling results of the third and forth case studies exhibit that evaporitic salinity conditions facilitate the development of full 3-D hydrothermal fluid flow system with the given physical parameters. Without the presence of the secondary cross fault, 3-D convection rolls tend to be stretched along the longitudinal strike direction of the major faults, and elongated discharge zones develop on the top of each major fault. With the presence of the cross fault, 3-D convection rolls tend to be ‘more’ rounded and characterised by ‘mushroom-like’ shapes. It is also clear from the comparison between Figs. 4c and 5b that the intersection of the two sets of the faults allows basinal fluids to channel downwards to a greater depth and circulate over larger volumes of the host rocks, and the heated basinal brines ascend to the basin floor via a more ‘localized’ discharge zone with greater temperature and venting velocity, which seems more favourable for the formation of giant ore deposits. 5. Discussions and conclusions 3-D numerical experiments have been conducted in this study in order to investigate the role of salinity variation in driving basin-scale hydrothermal ore-forming fluid flow. The general geological structure
Fig. 5. Scenario 4 corresponds to the condition that a cross fault is included to intersect the two major faults. Other conditions remain the same as in Scenario 3: (a) temperature contours at a time of 0.25 Ma; and (b) iso-surface of temperature distribution at the same time. Note that D and R denote the major discharge and recharge zones.
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
of the Mount Isa basin in northern Australia has been used to constrain our 3-D conceptualized model, that involves a sandstone aquifer and three other stratigraphical sequences cut by two steeply-dipping fault zones. Our numerical simulations have revealed that the interplay between the active synsedimentary faults and clastic aquifer unit controls the basinal fluid flow and hydrothermal discharge. Under open marine conditions in which salinity is assumed to remain constant throughout the basin, thermally-driven fluid flow circulates mainly within the more permeable faults themselves, rather than in the aquifer and other host rocks, unless the permeability contrast between the faults and aquifer is less than one order of magnitude. Our numerical experiments have also indicated that enhanced salinity on the basin floor resulting from evaporation facilitates buoyancy-driven free convection and can expel pore waters through discharge zone(s) to form SEDEX-type deposits at the basin floor by displacing lower salinity pore waters in the basin and starting 3-D thermohaline circulation of brines. If only the longitudinal major faults are present, 3-D hydrothermal convection rolls tend to be stretched along the fault strike direction so as to form elongated discharge zones on the basin floor. Including a cross fault to intersect the two primary major faults greatly alters basinal fluid flow patterns. The intersection of the two sets of the faults facilitates the development of ‘mushroom-shaped’ 3-D thermohaline convection rolls that allows the hot basinal brines to circulate through larger volumes of the host rocks to leach enough metal content, and more importantly it enables basinal fluids to channel downwards to a greater depth and eventually discharges to the basin floor via a more ‘localized’ discharge zone with higher temperature and venting velocity, which is favourable for the formation of giant SEDEX deposits. Geological implication of our numerical modeling results presented here is that SEDEX-type deposits are more easily formed when evaporation first produces surface brines and then these brines sink and displace pore waters in the basin. Also, determining the spatial distribution of fault intersection zones is critical to exploration for SEDEX deposits in sedimentary basins. It should be pointed out that this paper has not considered the role of tectonic deformation in driving ore-forming fluid migration; namely, presenting fluid flow patterns during the basin evolution of pre-structural inversion phase. Indeed, tectonic deformation and its effects on fluid flow are important for mineralization in northern Australia, in particular during the Isa orogenic inversion phase. Upon the commencement of tectonic deformation, fluid flow would switch from the convection patterns discussed here to more localized flow along the faults or in damage zones in the host rocks. Zhang et al. (2006) carried out the numerical modeling of fluid flow patterns during the shortening deformation event associated with the initial stages of the structural inversion in the Isa Superbasin at ca. 1575 Ma. Refer to their work for the details of deformation-driven fluid flow in the study area.
Acknowledgements This work was sponsored by the Program to Sponsor Teams for Innovation in the Construction of Talent Highlands in Guangxi Institutions of Higher Learning (during JY's sabbatical leave at the Guilin University of Technology, Guangxi, China). Research was also partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through a Discovery Grant to JY (Grant No. RGPIN 261283). The constructive comments from two anonymous reviewers are greatly appreciated.
243
References Barrier, C.T., Cathles, L.M., Erendi, A., 1999. Finite element heat and fluid flow computer simulations of a deep ultramafic sill model for the giant Kidd Creek volcanic associated massive sulfide deposit, Abititi Subprovince. Economic Geology Monograph 10, 529–540. Betts, P.G., Giles, D., Lister, G.S., 2003. Tectonic environment of shale-hosted massive sulfide Pb–Zn–Ag deposits of Proterozoic northeastern Australia. Economic Geology 98, 557–576. Bierlein, F.P., Betts, P.G., 2004. The Proterozoic Mount Isa Fault Zone, northeastern Australia: is it really a ca. 1.9 Ga terrane-bounding suture? Erath and Planetary Science Letters 225, 279–294. Feltrin, L., McLellan, J.G., Oliver, N.H.S., 2009. Modelling the giant, Zn–Pb–Ag Century deposit, Queensland, Australia. Computers & Geosciences 35, 108–133. 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–312. Hanor, J.S., 1996. Controls on the solubilization of lead and zinc in basin brines. In: Sangster, D.F. (Ed.), Carbonate-Hosted Lead–Zinc Deposits: Economic Geology Special Publication, vol. 4, pp. 483–500. Klein, S.A., Harvey, A.H., 1996. NIST/ASME Steam Properties (version 2.0), Standard Reference Data Program. National Institute of Standard and Technology, Maryland, USA. McLellan, J.G., Oliver, N.H.S., Hobbs, B.E., 2006. The relative effects of deformation and thermal advection on fluid pathways in basin-related mineralization. Journal of Geochemical Exploration 89, 271–275. Molson, J.W., Frind, E.O., Palmer, C.D., 1992. Thermal energy storage in an unconfined aquifer 2. model development, validation, and application. Water Resources Research 28, 2857–2867. O'dea, M.G., Lister, G.S., MacCready, T., 1997. Geodynamic evolution of the Proterozoic Mount Isa terrain. Geological Society of London Special Publication 121, 99–122. Oliver, N.H.S., McLellan, J.G., Hobbs, B.E., Cleverley, J.S., Ord, Al, Feltrin, L., 2006. Numerical models of extensional deformation, heat transfer, and fluid flow across basement– cover interfaces during basin-related mineralization. Economic Geology 101, 1–31. Sarkar, A., Nunn, J.A., Hanor, J.S., 1995. Free thermohaline convection beneath allochthonous salt sheets: an agent for salt dissolution and fluid flow in Gulf Coast sediments. Journal of Geophysical Research 100, 18,085–18,092. Schoofs, S., Spera, F.J., 2003. Transition to chaos and flow dynamics in thermochemical convection in porous media. Transport in Porous Media 50, 179–195. Yang, J., 2006a. Finite element modeling of transient saline hydrothermal fluids in multi-faulted sedimentary basins: implications for ore-forming processes. Canadian Journal of Earth Sciences 43, 1331–1340. Yang, J., 2006b. Full 3-D numerical simulation of hydrothermal fluid flow in faulted sedimentary basins: example of the McArthur Basin, Northern Australia. Journal of Geochemical Exploration 89, 440–444. Yang, J., Bull, S., Large, R., 2004. Numerical investigation of salinity in controlling oreforming fluid transport in sedimentary basins: example of the HYC deposit, northern Australia. Mineralium Deposita 39, 622–631. Yang, J., Large, R., Bull, S., Scott, D., 2006. Basin-scale numerical modelling to test the role of buoyancy driven fluid flow and heat transport in the formation of stratiform Zn–Pb–Ag deposits in the northern Mt. Isa basin. Economic Geology 101, 1275–1292. Zhang, Y., Sorjonen-Ward, P., Ord, A., Southgate, P.N., 2006. Fluid flow during deformation associated with structural closure of the Isa Superbasin at 1575 Ma in the central and northern Lawn Hill Platform, Northern Australia. Economic Geology 101, 1293–1312. Zhang, Y., Gartrell, A., Underschultz, J.R., Dewhurst, D.N., 2009. Numerical modelling of strain localisation and fluid flow during extensional fault reactivation: implications for hydrocarbon preservation. Journal of Structural Geology 31, 315–327. Zhao, C., Muhlhaus, H.B., Hobbs, B.E., 1997. Finite element analysis of steady-state natural convection problems in fluid-saturated porous media heated from below. International Journal for Numerical and Analytical Methods in Geomechanics 21, 863–881. Zhao, C., Hobbs, B.E., Muhlhaus, H.B., 1999. Theoretical and numerical analyses of convective instability in porous media with upward throughflow. International Journal for Numerical and Analytical Methods in Geomechanics 23, 629–646. Zhao, C., Hobbs, B.E., Muhlhaus, H.B., Ord, A., Lin, G., 2003. Finite element modeling of threedimensional steady-state convection and lead/zinc mineralization in fluid-saturated rocks. Journal of Computational Methods in Sciences and Engineering 3, 73–89. Zhao, C., Hobbs, B.E., Ord, A., Kuhn, M., Muhlhaus, H.B., Peng, S., 2006. Numerical simulation of double-diffusion driven convective flow and rock alteration in threedimensional fluid-saturated geological fault zones. Computer Methods in Applied Mechanics and Engineering 195, 2816–2840. Zhao, C., Hobbs, B.E., Ord, A., 2008a. Investigating dynamic mechanisms of geological phenomena using methodology of computational geosciences: an example of equaldistant mineralization in a fault. Science in China Series D: Earth Sciences 51, 947–954. Zhao, C., Hobbs, B.E., Ord, A., Hornby, P., Peng, S., 2008b. Morphological evolution of three-dimensional chemical dissolution front in fluid-saturated porous media: a numerical simulation approach. Geofluids 8, 113–127.
Contents lists available at ScienceDirect
Journal of Geochemical Exploration j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j g e o ex p
Three-dimensional numerical modeling of salinity variations in driving basin-scale ore-forming fluid flow: Example from Mount Isa Basin, northern Australia Jianwen Yang a,b,⁎, Zuohai Feng a, Xianrong Luo a, Yuanrong Chen a a b
Sabbatical leave at Faculty of Earth Sciences, Guilin University of Technology, Guilin, Guangxi, 541004, PR China Department of Earth and Environmental Sciences, University of Windsor, Windsor, Ontario, Canada N9B 3P4
a r t i c l e
i n f o
Article history: Received 27 April 2009 Accepted 16 December 2009 Available online 28 December 2009 Keywords: Hydrothermal fluid flow Finite element modeling SEDEX deposits Mount Isa basin
a b s t r a c t This paper presents a fully 3-D numerical investigation into the effect of salinity on ore-forming hydrothermal fluid flow, and develops a highly conceptualized hydrological model to simulate the fluid plumbing system early in the history of the Mount Isa basin, northern Australia when lead–zinc deposits were formed therein. Our numerical modeling results indicate that the active synsedimentary faults and clastic aquifer form a favourable hydro-framework for regional-scale fluid flow, and variations in salinity have important implications for fluid migration and heat transport. When salinity is constant throughout the basin, hydrothermal fluid flow mainly circulates within the more permeable faults as 2-D convection cells, unless the contrast in permeability between the faults and aquifer is less than one order of magnitude. Enhanced salinities on the basin floor due to evaporation of seawater facilitate the development of full 3-D thermohaline flow systems. The 3-D convection rolls established by the evaporitic conditions tend to be stretched if only longitudinal faults exist, but become ‘mushroom-shaped’ when a cross fault is added to intersect the longitudinal ones. The intersection of the two sets of the faults enables seawater to channel downwards to a greater depth and circulate through larger volumes of the source rocks to leach more metal content, and eventually the heated basinal brines ascend to the basin floor via a more ‘localized’ discharge zone at higher venting fluid velocity and temperature, potentially forming giant ore deposits. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Salinity (evaporitic) conditions of sedimentation have long been recognized as an important factor in the development of basinal brines that form stratiform Zn–Pb–Ag ore deposits, mainly because the solubility of the economically important metals increases with increasing salinity (e.g., Hanor, 1996). Also, spatial and temporal changes in salinity can lead to a buoyancy force due to density variation of basinal brine, which is critical in determining regional fluid flow patterns. Buoyancy-driven fluid flow has been extensively studied analytically and numerically in large-scale geological systems (e.g., Zhao et al., 1997, 1999). It is recently gaining momentum as the most likely hydrological scenario for the formation of sedimentary-exhalative (SEDEX) ore deposits. Previous numerical studies usually considered thermallyinduced fluid flow only, with no account taken of the effects of salinity (e.g., Barrier et al., 1999; Garven et al., 2001; Yang et al., 2006). A few studies have indeed considered both the heat- and salinity-induced buoyancy effects on fluid flow in porous media that are heated and ⁎ Corresponding author. Department of Earth and Environmental Sciences, University of Windsor, Windsor, Ontario, Canada N9B 3P4. Tel.: +1 519 253 3000x2181; fax: +1 519 973 7081. E-mail address: [email protected] (J. Yang). 0375-6742/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.gexplo.2009.12.004
salted from below due to the co-presence of natural geothermal gradients and salt domes (e.g., Sarkar et al., 1995; Schoofs and Spera, 2003), but they are not related to SEDEX ore genesis. We have recently developed a hydrogeological model that fully couples transient fluid flow, heat and solute transport associated with the formation of the HYC SEDEX deposit in the McArthur Basin, northern Australia (Yang et al., 2004; Yang, 2006a). However, these numerical experiments are 2-D in nature over the selected geological sections, and therefore cannot fully represent the complexities of 3-D hydrothermal system in reality. Recent fast developments in computer hardware and software have enabled researchers to start with simulating hydrothermal ore-forming fluid flow in realistic 3-D environments. For instance, Zhao et al. (2003) employed the finite element model to simulate 3-D steady-state convective pore-fluid flow and the relevant mineralization in fluidsaturated rocks. Zhao et al. (2006) derived analytical and numerical solutions for the double-diffusion driven convective instability in 3-D permeable faults, once again under steady-state conditions. Yang (2006b) presented a 3-D hydrological model that fully couples transient fluid flow and heat transport in the McArthur basin, northern Australia, but without considering the effect of salinity distribution on regional fluid flow. Zhao et al. (2008a) employed numerical modeling to address the controlling formation mechanism related to equal-distant distributed gold deposits in a 3-D fault. Zhao et al. (2008b) considered the morphological evolution of 3-D chemical dissolution fronts occurring in
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
fluid-saturated porous media. More recently Feltrin et al. (2009) simulated 3-D fluid flow driven by tectonic deformation associated with the formation of the giant lead–zinc–silver Century deposit, but no account was taken of the effect of temperature and salinity variation on fluid flow in the hydrothermal system. This paper presents a 3-D numerical investigation that fully couples transient fluid flow and heat transfer with solute transport, and it concerns in particular with buoyancy-driven fluid flow pertinent to the genesis of shale-hosted lead–zinc ores in the Mount Isa basin, northern Australia. In this study, we do not examine any specific deposits; rather we develop a highly generalized 3-D conceptual model on the basis of the general geological structure of the targeted basin. In addition in order to simplify the physics involved, we only consider a single solute transport and ignore chemical reactions between different phases. Ore-forming fluid flow can be driven by tectonic deformation as well (e.g., Zhang et al., 2009), however we do not intend to address this issue as our attention is focused on highlighting the role of salinity variations in driving basinal brines. We conduct a series of numerical experiments with different conditions in order to address the effect of salinity variations in driving hydrothermal ore-forming fluid migration in the study region. 2. Governing equations and finite element modeling Mathematically, saline hydrothermal fluid flow in sedimentary basins is governed by the fluid continuity equation, thermal energy conservation equation, and solute mass conservation equation, coupled through Darcy's law. The fluid continuity equation can be expressed in terms of an ‘equivalent freshwater’ hydraulic head for the variable density fluid flow system (Yang, 2006a). For a 3-D system in the x–y–z coordinate, we have: ∂ ∂h ∂ ∂h ∂ ∂h ∂h K + K + K + Kρr = Ss ; ∂x ∂x ∂y ∂y ∂z ∂z ∂t ∂ ∂T ∂ ∂T ∂ ∂T λm + λm + λm ∂x ∂x ∂y ∂y ∂z ∂z −
∂ ∂ ∂ ∂T ðc ρ q TÞ− ðcw ρw qx TÞ− ðcw ρw qz TÞ = cm ρm ; ∂x w w x ∂y ∂z ∂t
ð1Þ
ð2Þ
and ∂ ∂C ∂ ∂C ∂ ∂C θDx + θDy + θDz ∂x ∂x ∂y ∂y ∂z ∂z −
∂ ∂ ∂ ∂ ðq CÞ− ðqy CÞ− ðqz CÞ = ðθRCÞ; ∂x x ∂y ∂z ∂t
ð3Þ
where qx, qy and qz are the Darcy flux components in x, y, and z directions, h is the ‘equivalent freshwater’ head, K is the hydraulic − θ) conductivity, Ss is the specific storage, λm = λθwλ(1 (λw and λs s denote the thermal conductivity of the fluid and solid phase), θ is the porosity, cmρm = cwρwθ + csρs(1 − θ) (cw and cs are the specific heat capacity of the fluid and solid phase, and ρw and ρs are the density of the fluid and solid phase), ρr is the relative fluid density defined as ρr = ρw / ρ0 − 1 (ρ0 is the reference fresh water density), t is the time, C is the solute concentration, T is the temperature, R is the retardation factor, and Dx, Dy and Dz are the hydrodynamic dispersion coefficients. The governing Eqs. (1)–(3) are not alone sufficient. The dependence of fluid density and viscosity on temperature and solute concentration must be defined. In this paper it is calculated using the NIST/ASME Steam Properties code (Klein and Harvey, 1996). These equations form a time-dependent, nonlinear and coupled system, which leads their solutions to become nontrivial even for a very simple 1-D geological system. We have recently developed a finite element computer model, following original works by Molson et al. (1992), to numerically solve these governing equations. Listed
237
here are the major modeling steps: 1) update the fluid density and viscosity using the latest solute concentration and temperature from the previous time step; 2) determine the ‘equivalent freshwater’ head based on the updated fluid properties by solving the fluid continuity equation; 3) determine the Darcy flux based on Darcy's law using the updated fluid properties and hydraulic head; 4) calculate the temperature and solute concentration by respectively solving the thermal energy and solute mass conservation equations; and 5) check whether the preset convergence criteria are satisfactory. If not, steps 1) through 4) are repeated; if yes, move on to the next time step. We adopt a non-orthogonal quadrilateral mesh to create finite elements since it is better suited for the complex geometry of different stratigraphic units and faults encountered in sedimentary basins. Each element is assigned permeability, porosity, thermal conductivity and other physical parameters governing fluid flow, heat transport and solute transport based on the rock properties. Further details of the computational method, both the principles and the schemes for their implementation into finite element software, can be found from our previous publications (Yang et al., 2004; Yang, 2006a). 3. Three-dimensional conceptualized model As illustrated in Fig. 1a, our 3-D conceptual model has dimensions of 60, 30, and 20 km respectively in the x, y, and z directions, and is discretised by a 3-D non-orthogonal quadrilateral mesh consisting of 21,736 finite elements in total. Fig. 1b illustrates the paleohydrostratigraphy over the central x–z cross-section, which is constrained by some of the common features of a sedimentary basin's rift-and-sag phase, and in particular by the reconstructions of the Mount Isa basin, northern Australia (O'dea et al., 1997; Betts et al., 2003). Fig. 1 represents the highly conceptualized and simplified subsurface stratigraphy and structure, which controlled the hydrological system when the Mount Isa SEDEX deposits were formed early in the history of the Mount Isa basin. As shown in Fig. 1b, the conceptual model involves a volcanic basement sequence of low permeability (Unit 1), a sandstone aquifer of high permeability (Unit 2), a rift cover sequence of intermediate permeability (Unit 3), and an upper cover sequence of shales and siltstones (Unit 4) that hosts mineral deposits formed during or soon after sediment deposition. Two more permeable faults (Fault 1 and Fault 2) are also included to penetrate from the upper sequence into the basement. The faults are assumed to be 1 km wide and steeply-dipping, which is constrained by surface exposure and seismic profiling (Bierlein and Betts, 2004). A similar general geological structure of the Mount Isa–McArthur basin region has been also used in previous fluid flow models (e.g., Yang et al., 2004; McLellan et al., 2006; Oliver et al., 2006; Yang, 2006a). Like in the previous studies of McLellan et al. (2006) and Oliver et al. (2006), Fault 1 on the left is equivalent to the Mount Isa fault system, and Fault 2 on the right simulates a fault zone to farther north, such as the Termite Range fault at the Century deposit. The two steeplydipping faults cut the sandstone aquifer (Unit 2), forming a favourable hydrological framework for regional-scale fluid flow in both longitudinal and transverse directions. The permeabilities and thermal conductivities assigned to the stratigraphic elements and faults are given in Table 1, following the previous numerical investigations in the Mount Isa–McArthur basin region listed above. We assume that the vertical permeability of the host rocks is two orders of magnitude less than the horizontal permeability due to the bedded and stratified nature of these sedimentary rocks. In addition, we assume cw = 4174 J/kg °C (specific heat of the fluid phase), cs = 800 J/kg °C (specific heat of the solid phase), λw = 0.5 W/m °C (thermal conductivity of the fluid phase), ρs = 2630 kg/m3 (density of the solid phase), ρ0 = 1000 kg/m3 (reference fresh water density), and θ = 10% (porosity). The upper boundary (i.e., the basin floor) is permeable to fluid flow and its temperature is fixed at 20 °C. The lower boundary temperature is maintained at 450 °C, justified by present-day heat
238
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
Fig. 1. 3-D conceptualized hydrological model based on the general geological structure of the Mount Isa basin, northern Australia: (a) finite element mesh, and (b) central crosssection. Refer to Table 1 for the hydrological and thermal properties assigned to the faults and rock units.
flow measurements in northern Australia (S. McLaren and M. Sandiford, pers. comm.). The bottom of the model is assumed impermeable since it lies within the deep volcanic basement. The side boundaries are assumed to be adiabatic to heat transfer and impermeable to fluid flow. The initial fluid velocity is set to zero over the whole solution domain, and the initial temperature is assumed to vary linearly with depth. Various salinity conditions are considered in our numerical experiments as follows.
4. Numerical simulations and results The first scenario corresponds to a submarine condition, as geological evidence (e.g., Yang et al., 2006) indicates that open marine conditions persisted during deposition of the upper Isa Superbasin. We assume that basinal brine in this case is simply pure seawater, and therefore assign a constant salinity of typical seawater Table 1 Major physical parameters of the faults and host rocks. Hydrological unit Horizontal permeability Vertical permeability Thermal conductivity and formation (m2) (m2) (J/m/s/°C) Fault Unit 1 Unit 2 Unit 3 Unit 4
4 × 10− 14 2 × 10− 16 4 × 10− 14 2 × 10− 15 1 × 10− 15
4 × 10− 14 2 × 10− 18 4 × 10− 16 2 × 10− 17 1 × 10− 17
2.0 3.0 3.0 2.5 2.0
3.57% throughout the basin. Clearly under this condition, fluid flow is driven by thermally-induced buoyancy force only. Fig. 2 shows the numerical results corresponding to the reference conditions specified above and given in Fig. 1 and Table 1. Fig. 2a and b illustrates the temperature evolution from 0.5 to 0.8 Ma, and Fig. 2c shows the fluid velocity vectors at the time of 0.8 Ma. It can be seen that due to the large contrast of hydraulic conductivities between the faults and host rocks, significant fluid circulation is confined mainly within the more permeable faults with a maximum rate of 2.3 m/year rather than in the host media (Fig. 2c). Both upwelling and downwelling flow develops within these two faults at different longitudinal (y-axis) distances, which gives rise to significant temperature variation along the fault strike direction (Fig. 2b). By comparison, fluid flow in the host rocks including the aquifer is less significant, evidenced by the nearly flat isothermal lines (Fig. 2a and b). This is very similar to our previous 3-D modeling results of the McArthur basin under pure water conditions (Yang, 2006b). Since hydrothermal fluid flow more easily focuses and circulates within the fault zones with less fluid transport in the host rocks, the metal contents leached from the host rocks by the circulating basinal brines are likely insufficient to form giant SEDEX deposits. However, this situation may change if the sandstone aquifer (Unit 2) is more permeable. In the second scenario, therefore, we assume that the vertical permeability of the aquifer is increased by one order of magnitude (but still one order of magnitude lower than the permeability of the faults). Fig. 3 shows the temperature distribution and fluid velocity vectors at the time of 0.8 Ma, from which it can be
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
239
Fig. 2. Scenario 1 corresponds to a constant salinity of typical seawater of 3.57% assigned throughout the basin: (a) temperature contours at a time of 0.5 Ma; (b) temperature contours at a time of 0.8 Ma; and (c) fluid velocity vectors over the faults at a time of 0.8 Ma.
seen that two 3-D convection rolls now develop. Cold seawater recharges the basin via the central segment of Fault 1 (Fig. 3b and c), and it then not only circulates longitudinally along the fault like in the first scenario but also travels laterally through the aquifer toward the right (Fig. 3a and c). In the meantime, the basinal fluid is heated up from below and eventually ascends via the central segment of Fault 2 and discharges to the basin floor to form a SEDEX-type deposit (Fig. 3b and c). Fluid discharge temperatures at the surface of the model are in the range of 135 °C to 190 °C with fluid velocities of 1.5 to 2.5 m/year
for the period of 0.8 Ma. Note that as a 3-D convection roll, cold seawater also recharges the basin via the off-central parts of Fault 2 (Fig. 3a and b). It circulates downwards via the fault and the surrounding volume of the host rocks, heated up at depth, and finally flows upwards and discharges to the basin floor via the central part of Fault 2 (Fig. 3b and c). Our further numerical experiments (results not shown here) indicate that when the permeability contrast between the faults and aquifer is less than one order of magnitude, 3-D convective fluid flow system can always develop to circulate over a
240
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
Fig. 3. Scenario 2 corresponds to a constant salinity of typical seawater of 3.57% assigned throughout the basin, and the vertical permeability of the aquifer Unit 2 is increased by one order of magnitude from the reference value in Table 1: (a) temperature contours at a time of 0.8 Ma; (b) iso-surface of temperature distribution at a time of 0.8 Ma; and (c) fluid velocity vectors over the central cross-section at the same time. Note that D and R denote the major discharge and recharge zones.
large volume of the host rocks for leaching metals from the source rocks, which is in favour of the ore genesis of SEDEX-type deposits. However, when the permeability contrast exceeds one order of magnitude, fluid flow tends to be confined within the more permeable fault zone as 2-D convection cells, and this is not favourable for the ore formation. The third scenario corresponds to semi-evaporitic conditions resulting from the mixing of minimally modified seawater in the
sub-basins with saline brines generated by evaporation at the basin margins. We double the salinity of typical seawater 3.57% for the upper boundary (basin floor), and assign the seawater salinity for the lower boundary. Initial salinity is also equal to 3.57% throughout the basin. Fig. 4a and b shows the snapshots of iso-surface of the salinities 6% and 7% at times of 0.05 Ma and 0.25 Ma, corresponding to the reference conditions given in Fig. 1 and Table 1. It can be seen that this type of saline conditions facilitates buoyancy-driven free convection,
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
241
Fig. 4. Scenario 3 corresponds to the semi-evaporitic condition: (a) iso-surface of salinity distribution at a time of 0.05 Ma; (b) iso-surface of salinity distribution at a time of 0.25 Ma; and (c) iso-surface of temperature distribution at a time of 0.25 Ma. Note that D and R denote the major discharge and recharge zones.
compared with the pure seawater case discussed above. Due to the continuous supply of brine from the seafloor and hence the enhanced convection strength, saline fluid circulates not only within the faults but also spreads out via the aquifer and host rocks, and as a result the salinity in the aquifer increases progressively with time. Consequently several 3-D convection rolls develop, as identified from the iso-surface of temperatures 80 °C and 240 °C at the time of 0.25 Ma shown in Fig. 4c, which allows the hot basinal brines to circulate through sufficient volumes of the host rocks to leach enough metal content.
Also, as illustrated in Fig. 4c, several elongated discharge zones develop on the top of the faults, over which SEDEX-type deposits may be formed. The venting fluid temperatures on the basin floor range from 180 °C to 210 °C with fluid velocities of 47 to 75 m/year over the period of 0.25 Ma, both greater than those shown in the above scenarios under the pure seawater condition. In the forth scenario, we include a vertical cross fault to reflect the overall effect of several minor faults that sinuously and intermittently cut the longitudinal major faults in the targeted region (e.g., refer to
242
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
Fig. 1 of Yang et al., 2006). Other conditions remain the same as in the third scenario. The cross fault, located in the central x–z plane at y = 15 km, orthogonally intersects the two major faults. To simplify the complexity, this cross fault is assumed to have the same width and hydraulic conductivity as those of the primary main faults. The intersection of these two sets of faults may result in localized zones of enhanced permeability; therefore, two subvertical columns, centered by the left-hand and right-hand intersections and having a crosssectional area of 3 km by 3 km, are included and assigned the same permeability of the faults. Fig. 5 illustrates the temperature distribution at the time of 0.25 Ma, from which it can be seen that fluids are mainly channelled downwards via the left-hand localized column and upwards through the right-hand conduit. This leads to the development of two major 3-D ‘mushroom-shaped’ convection rolls that are centered by these two columns of enhanced permeability, and two ‘minor’ 3-D convection rolls that are located close to the ends of the cross faults (Fig. 5b). Significant hydrothermal convection now takes place not only along the longitudinal direction as shown in Fig. 2 of the first case study, but also along the transverse direction through both the cross fault and the host rocks, as evidenced by several upwelling and downwelling thermal plumes distributed alternately over the cross fault (Fig. 5a). The F1-associated left intersection zone behaves as the major recharge pathway, whereas the F2-related right localized zone serves as the major discharge conduit. Cold seawater moves downwards along the left subvertical column, then travels laterally to the right through the cross fault and the aquifer (the host rock Unit 2)
and in the meantime it is heated from below, and eventually the hot basinal brines ascend via the right subvertical column and discharge to the basin floor to form a SEDEX-type deposit. The venting fluid temperatures on the basin floor range from 190 °C to 235 °C with fluid velocities of 66 to 138 m/year over the period of 0.25 Ma. Numerical modeling results of the third and forth case studies exhibit that evaporitic salinity conditions facilitate the development of full 3-D hydrothermal fluid flow system with the given physical parameters. Without the presence of the secondary cross fault, 3-D convection rolls tend to be stretched along the longitudinal strike direction of the major faults, and elongated discharge zones develop on the top of each major fault. With the presence of the cross fault, 3-D convection rolls tend to be ‘more’ rounded and characterised by ‘mushroom-like’ shapes. It is also clear from the comparison between Figs. 4c and 5b that the intersection of the two sets of the faults allows basinal fluids to channel downwards to a greater depth and circulate over larger volumes of the host rocks, and the heated basinal brines ascend to the basin floor via a more ‘localized’ discharge zone with greater temperature and venting velocity, which seems more favourable for the formation of giant ore deposits. 5. Discussions and conclusions 3-D numerical experiments have been conducted in this study in order to investigate the role of salinity variation in driving basin-scale hydrothermal ore-forming fluid flow. The general geological structure
Fig. 5. Scenario 4 corresponds to the condition that a cross fault is included to intersect the two major faults. Other conditions remain the same as in Scenario 3: (a) temperature contours at a time of 0.25 Ma; and (b) iso-surface of temperature distribution at the same time. Note that D and R denote the major discharge and recharge zones.
J. Yang et al. / Journal of Geochemical Exploration 106 (2010) 236–243
of the Mount Isa basin in northern Australia has been used to constrain our 3-D conceptualized model, that involves a sandstone aquifer and three other stratigraphical sequences cut by two steeply-dipping fault zones. Our numerical simulations have revealed that the interplay between the active synsedimentary faults and clastic aquifer unit controls the basinal fluid flow and hydrothermal discharge. Under open marine conditions in which salinity is assumed to remain constant throughout the basin, thermally-driven fluid flow circulates mainly within the more permeable faults themselves, rather than in the aquifer and other host rocks, unless the permeability contrast between the faults and aquifer is less than one order of magnitude. Our numerical experiments have also indicated that enhanced salinity on the basin floor resulting from evaporation facilitates buoyancy-driven free convection and can expel pore waters through discharge zone(s) to form SEDEX-type deposits at the basin floor by displacing lower salinity pore waters in the basin and starting 3-D thermohaline circulation of brines. If only the longitudinal major faults are present, 3-D hydrothermal convection rolls tend to be stretched along the fault strike direction so as to form elongated discharge zones on the basin floor. Including a cross fault to intersect the two primary major faults greatly alters basinal fluid flow patterns. The intersection of the two sets of the faults facilitates the development of ‘mushroom-shaped’ 3-D thermohaline convection rolls that allows the hot basinal brines to circulate through larger volumes of the host rocks to leach enough metal content, and more importantly it enables basinal fluids to channel downwards to a greater depth and eventually discharges to the basin floor via a more ‘localized’ discharge zone with higher temperature and venting velocity, which is favourable for the formation of giant SEDEX deposits. Geological implication of our numerical modeling results presented here is that SEDEX-type deposits are more easily formed when evaporation first produces surface brines and then these brines sink and displace pore waters in the basin. Also, determining the spatial distribution of fault intersection zones is critical to exploration for SEDEX deposits in sedimentary basins. It should be pointed out that this paper has not considered the role of tectonic deformation in driving ore-forming fluid migration; namely, presenting fluid flow patterns during the basin evolution of pre-structural inversion phase. Indeed, tectonic deformation and its effects on fluid flow are important for mineralization in northern Australia, in particular during the Isa orogenic inversion phase. Upon the commencement of tectonic deformation, fluid flow would switch from the convection patterns discussed here to more localized flow along the faults or in damage zones in the host rocks. Zhang et al. (2006) carried out the numerical modeling of fluid flow patterns during the shortening deformation event associated with the initial stages of the structural inversion in the Isa Superbasin at ca. 1575 Ma. Refer to their work for the details of deformation-driven fluid flow in the study area.
Acknowledgements This work was sponsored by the Program to Sponsor Teams for Innovation in the Construction of Talent Highlands in Guangxi Institutions of Higher Learning (during JY's sabbatical leave at the Guilin University of Technology, Guangxi, China). Research was also partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through a Discovery Grant to JY (Grant No. RGPIN 261283). The constructive comments from two anonymous reviewers are greatly appreciated.
243
References Barrier, C.T., Cathles, L.M., Erendi, A., 1999. Finite element heat and fluid flow computer simulations of a deep ultramafic sill model for the giant Kidd Creek volcanic associated massive sulfide deposit, Abititi Subprovince. Economic Geology Monograph 10, 529–540. Betts, P.G., Giles, D., Lister, G.S., 2003. Tectonic environment of shale-hosted massive sulfide Pb–Zn–Ag deposits of Proterozoic northeastern Australia. Economic Geology 98, 557–576. Bierlein, F.P., Betts, P.G., 2004. The Proterozoic Mount Isa Fault Zone, northeastern Australia: is it really a ca. 1.9 Ga terrane-bounding suture? Erath and Planetary Science Letters 225, 279–294. Feltrin, L., McLellan, J.G., Oliver, N.H.S., 2009. Modelling the giant, Zn–Pb–Ag Century deposit, Queensland, Australia. Computers & Geosciences 35, 108–133. 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–312. Hanor, J.S., 1996. Controls on the solubilization of lead and zinc in basin brines. In: Sangster, D.F. (Ed.), Carbonate-Hosted Lead–Zinc Deposits: Economic Geology Special Publication, vol. 4, pp. 483–500. Klein, S.A., Harvey, A.H., 1996. NIST/ASME Steam Properties (version 2.0), Standard Reference Data Program. National Institute of Standard and Technology, Maryland, USA. McLellan, J.G., Oliver, N.H.S., Hobbs, B.E., 2006. The relative effects of deformation and thermal advection on fluid pathways in basin-related mineralization. Journal of Geochemical Exploration 89, 271–275. Molson, J.W., Frind, E.O., Palmer, C.D., 1992. Thermal energy storage in an unconfined aquifer 2. model development, validation, and application. Water Resources Research 28, 2857–2867. O'dea, M.G., Lister, G.S., MacCready, T., 1997. Geodynamic evolution of the Proterozoic Mount Isa terrain. Geological Society of London Special Publication 121, 99–122. Oliver, N.H.S., McLellan, J.G., Hobbs, B.E., Cleverley, J.S., Ord, Al, Feltrin, L., 2006. Numerical models of extensional deformation, heat transfer, and fluid flow across basement– cover interfaces during basin-related mineralization. Economic Geology 101, 1–31. Sarkar, A., Nunn, J.A., Hanor, J.S., 1995. Free thermohaline convection beneath allochthonous salt sheets: an agent for salt dissolution and fluid flow in Gulf Coast sediments. Journal of Geophysical Research 100, 18,085–18,092. Schoofs, S., Spera, F.J., 2003. Transition to chaos and flow dynamics in thermochemical convection in porous media. Transport in Porous Media 50, 179–195. Yang, J., 2006a. Finite element modeling of transient saline hydrothermal fluids in multi-faulted sedimentary basins: implications for ore-forming processes. Canadian Journal of Earth Sciences 43, 1331–1340. Yang, J., 2006b. Full 3-D numerical simulation of hydrothermal fluid flow in faulted sedimentary basins: example of the McArthur Basin, Northern Australia. Journal of Geochemical Exploration 89, 440–444. Yang, J., Bull, S., Large, R., 2004. Numerical investigation of salinity in controlling oreforming fluid transport in sedimentary basins: example of the HYC deposit, northern Australia. Mineralium Deposita 39, 622–631. Yang, J., Large, R., Bull, S., Scott, D., 2006. Basin-scale numerical modelling to test the role of buoyancy driven fluid flow and heat transport in the formation of stratiform Zn–Pb–Ag deposits in the northern Mt. Isa basin. Economic Geology 101, 1275–1292. Zhang, Y., Sorjonen-Ward, P., Ord, A., Southgate, P.N., 2006. Fluid flow during deformation associated with structural closure of the Isa Superbasin at 1575 Ma in the central and northern Lawn Hill Platform, Northern Australia. Economic Geology 101, 1293–1312. Zhang, Y., Gartrell, A., Underschultz, J.R., Dewhurst, D.N., 2009. Numerical modelling of strain localisation and fluid flow during extensional fault reactivation: implications for hydrocarbon preservation. Journal of Structural Geology 31, 315–327. Zhao, C., Muhlhaus, H.B., Hobbs, B.E., 1997. Finite element analysis of steady-state natural convection problems in fluid-saturated porous media heated from below. International Journal for Numerical and Analytical Methods in Geomechanics 21, 863–881. Zhao, C., Hobbs, B.E., Muhlhaus, H.B., 1999. Theoretical and numerical analyses of convective instability in porous media with upward throughflow. International Journal for Numerical and Analytical Methods in Geomechanics 23, 629–646. Zhao, C., Hobbs, B.E., Muhlhaus, H.B., Ord, A., Lin, G., 2003. Finite element modeling of threedimensional steady-state convection and lead/zinc mineralization in fluid-saturated rocks. Journal of Computational Methods in Sciences and Engineering 3, 73–89. Zhao, C., Hobbs, B.E., Ord, A., Kuhn, M., Muhlhaus, H.B., Peng, S., 2006. Numerical simulation of double-diffusion driven convective flow and rock alteration in threedimensional fluid-saturated geological fault zones. Computer Methods in Applied Mechanics and Engineering 195, 2816–2840. Zhao, C., Hobbs, B.E., Ord, A., 2008a. Investigating dynamic mechanisms of geological phenomena using methodology of computational geosciences: an example of equaldistant mineralization in a fault. Science in China Series D: Earth Sciences 51, 947–954. Zhao, C., Hobbs, B.E., Ord, A., Hornby, P., Peng, S., 2008b. Morphological evolution of three-dimensional chemical dissolution front in fluid-saturated porous media: a numerical simulation approach. Geofluids 8, 113–127.