Some fundamental issues in computational hydrodynamics of mineralization: A review

Journal of Geochemical Exploration 112 (2012) 21–34

Contents lists available at SciVerse ScienceDirect

Journal of Geochemical Exploration journal homepage: www.elsevier.com/locate/jgeoexp

Some fundamental issues in computational hydrodynamics of mineralization: A review Chongbin Zhao a,⁎, Lynn B. Reid b, c, Klaus Regenauer-Lieb b, d a

Computational Geosciences Research Centre, Central South University, Changsha 410083, China Western Australian Geothermal Centre of Excellence, CSIRO Division of Earth Science and Resource Engineering, P. O. Box 1130, Bentley, WA 6102, Australia School of Environmental System Engineering, The University of Western Australia, Crawley, WA 6009, Australia d School of Earth and Environment, The University of Western Australia, Crawley, WA 6009, Australia b c

a r t i c l e

i n f o

Article history: Received 18 October 2010 Accepted 28 October 2011 Available online 4 November 2011 Keywords: Fundamental issue Fluid flow Ore forming system Hydrodynamics Mineralization Computational simulation

a b s t r a c t This paper presents a general discussion on some fundamental issues associated with the current research status in the field of computational hydrodynamics of mineralization. Based on the scientific characteristics of mineral forming systems within the upper crust of the Earth, the fundamental issues under consideration include: (1) simulation of multi-process aspects of a mineral forming system; (2) simulation of multi-scale aspects of a mineral forming system; (3) simulation of fluid convection associated with a mineral forming system; (4) simulation of fluid focusing and mixing associated with a mineral forming system; and (5) simulation of reactive mass transport associated with a mineral forming system. After the current status of each fundamental issue is briefly summarized, the future development in the field of computational hydrodynamics of mineralization is discussed on the basis of further improving the treatments of the fundamental issues considered in this investigation. © 2011 Elsevier B.V. All rights reserved.

Contents 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2. Simulation of multi-process aspects of a mineral forming system . . . . . . . . . 3. Simulation of multi-scale aspects of a mineral forming system . . . . . . . . . 4. Simulation of fluid convection associated with a mineral forming system . . . . 5. Simulation of fluid focusing and mixing associated with a mineral forming system 6. Simulation of reactive mass transport associated with a mineral forming system . 7. Conclusions and discussions . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1. Introduction With the ever-increasing demand for improving our living standards, it has been recognized that existing natural resources will be exhausted in the near future and that our living environments are, in fact, deteriorating as a result of the greenhouse effect. With mining industries in Australia and China taken as an example, most of the

⁎ Corresponding author. E-mail address: [email protected] (C. Zhao). 0375-6742/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.gexplo.2011.10.005

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

21 23 24 27 28 29 31 32 32

easy-to-mine, near-surface mineral deposits have been discovered and developed. Although under-explored regions such as west Africa, eastern Europe, Latin America and south-east Asia still have lots of untapped potential, fundamental mining infrastructures are generally lacking in these regions. This means that huge financial investments are needed to explore and extract mineral resources in such underexplored regions. To make best use of existing mining infrastructures so as to reduce mining costs significantly, it is reasonable to explore concealed resources in the surrounding areas of existing mines. Since a large amount of geological data is already available as a result of an existing mine excavation, it is possible to investigate the major physical and chemical processes that control the formation of ore

22

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

deposits in the region where an existing mine is located. This means that in order to maintain the sustainable development of our living standards and the further improvement of our living environments, an inevitable and challenging task that geoscientists are now confronting is how accurately to predict not only the occurrences of natural disasters such as earthquakes, but also the locations of large concealed natural resources in the deep Earth. For this reason, geoscientists must study the processes, rules and laws, by which the Earth system operates, instead of simply describing and observing geoscience phenomena. On the other hand, with tremendous advances in contemporary computer technology during the past few decades, numerical methods and computational simulations provide a new way to deal with many geoscience problems. Traditional theoretical and experimental methods may not be valid in some problems as a result of the large time and length scales involved. Many hitherto unsolvable geoscience problems have been solved using numerical methods and computational simulations (Garven and Freeze, 1984; Raffensperger and Garven, 1995; Schafer et al., 1998; Steefel and Lasaga, 1994; Yeh and Tripathi, 1991; Zhao et al., 1997, 2008a, 2009). In particular, through wide application of computational science to geoscience problems, a new discipline, namely computational geoscience, has been established (Zhao et al., 2009). Since computational geoscience is established through applying a combination of mathematics, physics, chemistry and advanced computer algorithms to solve geoscience problems, it is obviously of multi-disciplinary nature across many fields of science. For this reason, the methodology of computational geoscience is a comprehensive research methodology, which is formed by integrating field observation, theoretical analysis, computational simulation and field validation (Zhao et al., 2008b, 2009). The primary aim of using the methodology of computational geoscience is to investigate dynamic processes and mechanisms of observed geological phenomena, rather than to describe the observed geological phenomena themselves. Although the emerging computational geoscience discipline can include many subdisciplines, such as computational geodynamics, computational geology, computational geochemistry, computational geomechanics, computational geophysics, computational hydrogeology and so forth, the main focus of this paper is on the current developing status and several fundamental issues associated with computational hydrodynamics of mineralization. Since the main purpose of computational hydrodynamics of mineralization is to simulate ore forming systems within the upper crust of the Earth, the major physical and chemical processes to be considered should include the following controlling processes: the material deformation process, fluid flow process, heat transfer process, mass transport process and chemical reaction process (Garven and Freeze, 1984; Raffensperger and Garven, 1995; Schafer et al., 1998; Schaubs and Zhao, 2002; Steefel and Lasaga, 1994; Yeh and Tripathi, 1991; Zhao et al., 1997, 1998, 2006, 2007a). As these processes are often fully coupled in a nonlinear, multiscale manner, we must not only develop modeling capabilities to simulate each of the above processes, but also, we have to develop the ability to simulate the forward and backward interactions among them in a seamless, coordinated manner. Using the methodology of computational geoscience, the following three kinds of models need to be established for a typical mineralization problem: a conceptual model, a mathematical model and a computational simulation model. In terms of establishing a conceptual model, the following five questions need to be answered (Gow et al., 2002; Hobbs et al., 2000; Ord et al., 2002; Schaubs and Zhao, 2002; Sorjonen-Ward et al., 2002): (1) What is the geometry and material distribution characteristics of the entire mineralizing system? (2) What is the geological history of the system with an emphasis on timing, pressure and temperature history? (3) What processes were responsible for driving fluid flow? (4) What were the physical and chemical processes involved with (and signatures of) the fluids responsible for mineralization? (5) What were the processes involved in the precipitation and development of minerals and

alteration at the site of mineralization? The answers to these five questions are fundamental to the construction of realistic conceptual models of a mineralizing system and to supply detail of geometrical size, material distributions, physical and chemical processes, initial and boundary conditions for these models. In terms of establishing a mathematical model, three fundamental scientific principles, namely the conservation laws of mass, momentum and energy, as well as the related physical and chemical theorems are commonly used to describe the controlling processes associated with the mineralizing system within the upper crust of the Earth. As a result, a set of partial differential equations is obtained from the derivation of the mathematical model. Using numerical discretization methods such as the finite element or finite difference method, a computational simulation model consisting of a set of algebraic equations can be established from the derived mathematical model. After this set of algebraic equations is solved, numerical solutions of an approximate nature can be obtained for the considered mineralizing system. Since numerical methods produce approximate solutions for ore forming systems, computational simulation models work well only if they employ appropriate numerical methods and algorithms, and their limitations are appropriately appreciated. For example, the numerical solution reliability of a mineralizing system is strongly dependent on the algorithm convergence, algorithm stability, mesh shape, time-step and other factors. To ensure the accuracy and reliability of computational simulation results for an ore forming system, the above-mentioned factors need to be carefully considered in the process of establishing computational simulation models. A newly-developed computer program needs to be verified through corresponding benchmark problems (Zhao et al., 1997) before it is used to simulate any real mineralizing systems. Otherwise, the reliability of the computational simulation solution obtained from the newly-developed computer program cannot be guaranteed. Since a mineralizing system within the upper crust of the Earth is commonly regarded as a complex system with multi-scale and multiprocess characteristics, it is still very difficult to simultaneously simulate the entire mineralizing system, from a particle (i.e. mineral) scale to a crustal scale, in one computational model. Nevertheless, it is possible to separately simulate different key aspects of a mineralizing system using a series of computational models consisting of different processes and different scales. This compartmentalization allows only some key factors to be considered for each of the computational models. By integrating all the simulation results obtained from a series of computational models consisting of varying processes and scales, it is possible to understand the formation mechanisms behind mineralizing systems. From the continuum mechanics point of view, the crustal rock can be treated as a porous medium, in which porosity is used to describe the pore volume per unit volume rock. In such a porous medium, fluid can flow in the flow channels that are formed by the connected pores and fractures. On the other hand, volumetric strain is used to describe the volume change per unit volume rock as a result of rock deformation. Since solid particles in a porous medium are usually assumed to be rigid (i.e. uncompressible), a porosity change is related to the corresponding volumetric strain change during a deformation process. Thus, any volumetric strain change can result in a change in the microscopic flow channels through the corresponding change in porosity. Due to the complex and complicated nature of a mineralizing system, only some fundamental issues of common nature in the computational simulation of mineralizing systems are addressed in this paper. Toward this end, other fundamental issues, such as multiphase flow simulation, magma intrusion process simulation and parameters selection, are not considered, even though they are important for the computational simulation of several specific mineralizing systems. It must be pointed out that although multiphase flow simulation has been well developed and used in other research fields (Bergins et al., 2005; Crone et al., 2002; Driesner and Geiger, 2007; Miller et al.,

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

1998; Schulenberg and Muller, 1987; Wang and Beckermann, 1993; Wang and Cheng, 1996), it is limited in the computational simulation of mineralizing systems (Geiger et al., 2006a,b; Mills et al., 2007, 2009). Similarly, there is little work available in the computational simulation of magma intrusion processes, although the thermal and chemical effects of magma intrusion on the fluid flow and mineralization have been investigated (Zhao et al., 2009). On the other hand, because literature materials available are extensive for the fundamental issues under consideration, we have had to be highly selective and have chosen references that address the specific points of an issue rather than attempt to be inclusive. Thus, the main purpose of this paper is to address some fundamental issues of common nature in the computational simulation of mineral forming systems, rather than to review the existing literature that is associated with the computational simulation of various types of mineral deposits. For this purpose, only generic mineralizing systems consisting of different processes and different scales are used, in this investigation, to discuss some fundamental issues in the field of computational hydrodynamics of mineralization. The fundamental issues to be considered in this paper are: (1) simulation of multi-process aspects of a mineral forming system; (2) simulation of multi-scale aspects of a mineral forming system; (3) simulation of fluid convection associated with a mineral forming system; (4) simulation of fluid focusing and mixing associated with a mineral forming system; and (5) simulation of reactive mass transport associated with a mineral forming system. 2. Simulation of multi-process aspects of a mineral forming system The central topic of the computational hydrodynamics of mineralization is to consider the contributions of fluid flow to ore body formation within the upper crust of the Earth. Since crustal materials are usually comprised of porous rocks, not only can mineral bearing fluids reside in the pore space of a porous rock, but also they can flow through the connected pores in the porous rock. Thus, fluid flow becomes the sole agent to transport dissolved-minerals from one place to another within the Earth's upper crust (Schmidt Mumm et al., 2010; Zhao et al., 2010). Generally, fluid flow can be generated by mechanical, thermal or chemical processes, or combinations of these. In the case of mechanical processes, fluid flow mainly results from a deformation or gravity force induced pressure-gradient and therefore is called pressure-gradient driven flow. This kind of fluid flow can cause advective heat transfer and mass transport in the fluid-saturated porous rock. Fluid flow induced by the uneven surface topography of a basin and the flow of fluids squeezed out of a sedimentary basin by crustal deformation are typical examples of fluid pressure-gradient driven flow (Bethke, 1985; Ge and Garven, 1992; Oliver, 1986). In the case of thermal processes, fluid flow is predominantly caused by a temperature gradient so that it is often called thermal convection. Similarly, this kind of fluid flow can lead to convective heat transfer and mass transport in fluid-saturated porous rocks (Nield and Bejan, 1992; Zhao et al., 1997, 1998). In the case of chemical processes, the density effect of a chemical-species concentration gradient is the main driving force to trigger fluid flow. This kind of fluid flow is often called concentration gradient driven flow (Nield and Bejan, 1992; Yang et al., 2010; Zhao et al., 2008a). Note that although a chemical-species concentration gradient can cause mass diffusion, which is a transport process but not a flow process, the fluid density variation arising from a chemical-species concentration gradient can cause fluid flow. To determine the controlling processes that need to be considered for mineralizing systems involved in the field of computational hydrodynamics of mineralization, it is necessary to identify the key factors that may affect both the fluid attributes (such as density, viscosity, conductivity and so forth) and the creation of flow channels in fluidsaturated porous rocks. Since a change in fluid pressure can cause a change in both the effective stress in a porous rock and the strength of the porous rock, it can cause a considerable change in the flow channels

23

of the porous rock through material deformation processes. This is because the effective stress in a porous rock is directly related to the volumetric strain through the constitutive law of the rock. Thus, any change in the effective stress of the rock can lead to a volumetric strain of the rock. As mentioned in the introduction, since any volumetric strain change can result in a change in the microscopic flow channels through the corresponding change in porosity, any change in the effective stress of the rock can cause a change in the microscopic flow channels of the porous rock through material deformation processes. On the other hand, fluid pressure can cause the initiation and propagation of cracks in porous rocks through hydraulic fracturing, which can be regarded as a particular kind of material deformation process. Note that hydraulic fracturing plays an important role in creating flow channels associated with magma intrusion problems in the crust of the Earth (Bons et al., 2001; Lister, 1990; Lister and Kerr, 1991; Marsh, 1982; Rubin, 1995; Scott and Stevenson, 1986; Weinberg, 1996; Zhao et al., 2008c). Although many material deformation phenomena resulting from tectonic activities such as earthquakes, faulting, thrusting and folding can create flow channels in the crust of the Earth, their time scales are much smaller than the time scale of mineralization, as discussed in the next section of this paper. Nevertheless, the geological structures created during tectonic activities can provide useful flow channels and therefore may play some controlling roles in ore body formation and mineralization during post-tectonic activities in the upper crust of the Earth. Corresponding to different contributions of material deformation to the computational hydrodynamics of mineralization, material deformation can be divided into the following two categories: premineralization deformation and syn-mineralization deformation. In the former case, material deformation (that can result in a change in the volumetric strain of the rock) takes place before mineralization and therefore provides initial flow channels for mineral-bearing fluids in a mineralizing system, while in the latter case, material deformation (that can result in a change in the volumetric strain of the rock) takes place just during mineralization and therefore can cause a change in flow channels for mineral-bearing fluids in a mineralizing system. Since most pre-mineralization deformation resulting from tectonic activities is often large and discontinuous, as in faulting, the particle simulation method has been proven useful for simulating such premineralization deformation problems occurring in the crust of the Earth (Burbidge and Braun, 2002; Donze et al., 1994; Finch et al., 2004; Owen et al., 2004; Saltzer and Pollard, 1992; Zhao et al., 2007b, 2008c). Fluid flow in fluid-saturated porous rocks is commonly described using Darcy's law, while magma flow in magma chambers is simulated using the Navier–Stokes equation. Generally, the density of fluid is dependent on temperature, pressure and the salinity (total dissolved species concentration) in the pore space, while the viscosity of fluid is mainly dependent on temperature (Lin et al., 2003; Zhao et al., 2008a). To facilitate mathematical treatment in a theoretical analysis, the equation of state of fluid is usually replaced by the Oberbeck– Boussinesq approximation (Nield and Bejan, 1992; Phillips, 1991; Zhao et al., 1997), in which the density of fluid is treated as either a function of temperature for thermal diffusion driven convection problems or a function of temperature and species concentration for double diffusion driven convection problems (Frind, 1982; Harcouet-Menou et al., 2009; Kuhn et al., 2006; Yang et al., 2010; Zhao et al., 2006). However, for the computational simulation of crustal-scale fluid flow, if possible, a realistic equation of state of fluid needs to be used to simultaneously consider the dependence of fluid density on temperature, pressure and salinity in the pore space. In terms of simulating heat transfer processes associated with a mineralizing system in the field of computational hydrodynamics of mineralization, Fourier's law is used to derive the energy balance equation of fluid in fluid-saturated porous rocks (Dagan, 1989; Nield and Bejan, 1992). Due to the slowness of fluid flow in the crust of the Earth, thermal dispersion, which is directly proportional to the velocity of fluid flow, is much smaller than thermal diffusion

24

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

in the fluid-saturated porous rock. This means that thermal dispersion can be safely neglected in the conceptual model. Also, the slowness of fluid flow enables viscous dissipation to be neglected in the computational simulation of heat transfer processes in mineralizing systems within the upper crust of the Earth. As a potential heat source in the crust of the Earth, radioactive heat may be generated by the decay of radioactive elements. Since many elements have unstable, radioactive isotopes, a small amount of heat is released when a parent isotope decays to its daughter isotope. However, only three radioactive elements (namely uranium, thorium and potassium), which occur within the Earth in sufficient abundance and release enough heat per decay event, are significant contributors to radioactive heat production (Clauser, 2009). Uranium has two radioactive isotopes (uranium-235 ( 235U) and uranium-238 ( 238U)); thorium and potassium have one each (thorium-232 ( 232Th) and potassium-40 ( 40 K)). Since each so-called heatproducing isotope decays at a different rate, the total rate of radioactive heat generation has declined and the relative importance of each isotope has varied over time. Uranium, thorium and potassium are all much more abundant in the crust than in the mantle, because they are all “incompatible” elements that go preferentially into the magma when rock begins to melt. The mantle-derived melts from which the Earth's crust formed are therefore greatly enriched in these elements. Nevertheless, radioactive heat generation is not uniform within the crust. It is the greatest in rocks of the granite family, in which the incompatible elements are most highly enriched. This implies that the rate of heat generation is the greatest in the upper part of the continental crust, where granites are the most abundant, and decreases with depth. As a result, the temperature gradient in the upper crust averages 30 °C/km (Clauser, 2009). Since an oceanic crust, which is basaltic in composition, is much poorer in uranium, thorium and potassium, less heat is produced there for an equal volume of rock. The generated heat in the oceanic crust is even less than that in the lower part of the continental crust. For the above reasons, when a mineralizing system is located in the continental crust, radioactive heat production (which is a transient process) needs to be considered in the computational simulation of the mineralizing system. Unfortunately, it is impossible to obtain an accurate estimate of the Earth's present-day rate of radioactive heat production because the total abundances of the heat-producing radioactive elements are poorly known. Although they can be measured precisely in any single sample taken, it is difficult to estimate a reliable global average because the crust is so variable in composition. Consequently, the rate of heat production even in the relatively welldocumented upper crust is uncertain by about 50% either way. For this reason, one may apply an appropriate geothermal gradient (as an initial condition) to approximately simulate the final effect of radioactive-heat generation in a computational model. This means that the transient process of radioactive heat production is neglected in the computational model. To examine how the transient process of radioactive heat production may affect ore forming processes, further research is needed to provide an appropriate radioactive heat production rate for the computational simulation of a mineralizing system in the upper part of a continental crust. Another potential heat source/sink relevant to ore body formation and mineralization comes from exothermic and endothermic chemical reactions involved in a mineralizing system. Since chemical reactions release or absorb energy, they can affect the temperature of their surrounding rocks. Exothermic reactions release energy to the surrounding rocks and cause them to heat up, while endothermic reactions absorb energy from the surrounding rocks and cause them to cool down. Due to the relative small amount of heat released from or absorbed by chemical reactions, this factor was generally neglected in conceptual models.

In order to simulate mass transport processes associated with a mineralizing system in the field of computational hydrodynamics of mineralization, Fick's law is used to derive the mass conservation equations of primary species in fluid-saturated porous rocks (Dagan, 1989). Generally, the total number of mass conservation equations under consideration is equal to the total number of primary species of interest in the fluid-saturated porous rocks. Unlike thermal dispersion, solute dispersion, which is directly proportional to the velocity of fluid flow, is usually more important than solute diffusion in fluid-saturated porous rocks so that it needs to be considered in the computational simulation of a mineralizing system (Dagan, 1989). Currently, there are two ways, namely an approximate way and an accurate way, to be used for considering solute dispersion in a computational model. In the approximate way, the solute dispersion effect is assumed to be of isotropic nature so that it can be roughly considered through increasing the solute diffusion coefficient in the fluid-saturated rock. To facilitate mathematical treatments, this technique has been commonly used in the theoretical analysis of mass transport problems in fluid-saturated porous rocks (Nield and Bejan, 1992; Phillips, 1991; Zhao et al., 2008a, 2010). In the accurate way, the solute dispersion effect is described using a dispersion tensor so that the anisotropic nature of solute dispersion can be simulated in a computational model. An important mass source/sink for a particular mineral under consideration is the dissolution or precipitation of this mineral through heterogeneous chemical reactions that take place in the considered mineralizing system. A homogeneous chemical reaction is referred to as the chemical reaction taking place within the same phase, such as between aqueous species in the fluid in the considered geochemical system, while a heterogeneous chemical reaction is referred to as the chemical reaction taking place between two or more phases, such as between the fluid and the minerals. As a typical example, a fluid–rock interaction problem usually involves at least one heterogeneous chemical reaction (Alt-Epping and Zhao, 2010; Zhao et al., 2001). To accurately determine the mass source/sink of a mineral, the thermodynamic properties of all chemical reactions involved in a mineralizing system under consideration have to be known. As quantitative modeling has become a standard procedure in geochemistry, there is a wealth of thermodynamic data available for a broad range of minerals, aqueous species and gases that allow the quantitative evaluation of mineral dissolution and precipitation (Barnes, 1997; Geiger et al., 2006a,b; Pokrovskii, 1999; Shock and Helgeson, 1988; Shock and Koretsky, 1995; Wang et al., 2010). 3. Simulation of multi-scale aspects of a mineral forming system As mentioned in the previous section, a mineralizing system to be considered in the field of computational hydrodynamics of mineralization may involve the following main processes: (1) material deformation process; (2) fluid flow process; (3) heat transfer process; (4) mass transport process and (5) chemical reaction process. To demonstrate the multiple time-scale characteristics of a mineralizing system, it is desirable to estimate the time-scale for each of the main processes involved in the computational simulation of the mineralizing system. When material deformation is caused by a mechanical process, the time-scale of the material deformation process can be estimated by considering stress diffusion in porous rocks. According to the classical elastic theory, the speed of the fastest propagating wave, namely the P-wave, in an elastic, homogeneous and isotropic medium can be expressed as follows:

VP ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi κ þ 43 G ρ

ð1Þ

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

where VP is the speed of the P-wave; κ is the bulk modulus of the medium; G is the shear modulus of the medium and ρ is the density of the medium. Note that both the bulk modulus and the shear modulus of an elastic medium can be also expressed in the following formulas: E ; κ¼ 3ð1−2ν Þ

E G¼ 2ð1 þ ν Þ

ð2Þ

where E is the elastic modulus (i.e. Yong's modulus) of the medium and ν is Poisson's ratio of the medium. Based on the P-wave speed of a stress wave in a rock mass, the characteristic time (i.e. the time-scale) of material deformation as a result of stress diffusion in the rock mass can be approximately estimated using the following expression: t deformation

Ldeformation ¼ ¼ Ldeformation VP

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ρð1−2ν Þð1 þ νÞ Eð1−ν Þ

ð3Þ

where tdeformation is the characteristic time of the material deformation process in a mineralizing system and Ldeformation is the characteristic length of the material deformation process, which may be equal to the crustal thickness of the mineralizing system. Since the P-wave speed is equal to a few thousands of kilometers per second in the Earth's crust and the thickness of a crust is a few tens of kilometers, the characteristic time of the material deformation process in a mineralizing system of the crustal length-scale may vary within a range of about a few milli-seconds. For example, if an upper crust has a thickness of 10 km and the P-wave speed is 4000 km/s, then the characteristic time of the considered system is equal to 0.0025 s. Note that to accurately simulate the material deformation process in the crust, the time-step used in a computational model should be much smaller than the characteristic time of the considered system. The time-scale of the fluid flow process involved in a mineralizing system can be estimated by considering fluid pressure diffusion in the fluid-saturated porous rocks. This can be achieved through investigating the specific formulation of a fluid pressure diffusion problem in the fluid-saturated porous rock. Toward this goal, it is necessary to derive the governing equation of fluid pressure diffusion in a fluidsaturated porous medium. If fluid is slightly compressible, the continuity equation of the fluid flow can be expressed in the Cartesian coordinate system as follows:

φ

∂ρf ∂t

þ ρf

  ∂u ∂v ∂w þ þ ¼0 ∂x ∂y ∂z

ð4Þ

where u, v and w are the Darcy velocities of fluid flow in the x, y and z directions respectively; ρf is the density of the fluid and φ is the porosity of the porous rock. According to Darcy's law, the related Darcy velocities are expressed in the following equations: u¼−

K x ∂p μ ∂x

ð5Þ

v¼−

K y ∂p μ ∂y

ð6Þ

K z ∂p μ ∂z

ð7Þ

w¼−

where p is the excess pressure of the fluid; μ is the dynamic viscosity of the fluid; Kx, Ky and Kz are the permeability of the porous medium in the x, y and z directions respectively. Note that the excess pressure of the fluid is equal to the total pressure of the fluid minus the

25

hydrostatic pressure of the fluid. If the porous medium is homogeneous and isotropic, then the following expression exists: K x ¼ K y ¼ K z ¼ K control

ð8Þ

where Kcontrol are the controlling permeability of the porous medium. In the case of an anisotropic porous medium, the controlling permeability can be determined by taking the maximum value among Kx, Ky and Kz. From the classical fluid dynamics (Gerhart et al., 1992), the bulk modulus of elasticity for a slightly compressible fluid can be expressed as follows: Ev ¼ −V f

dp dp ¼ ρf dV f dρf

ð9Þ

where Ev is the bulk modulus of the fluid; ρf is the density of the fluid and Vf is the volume of the fluid. Substituting Eqs. (5) to (9) into Eq. (4) yields the following diffusion equation of fluid pressure in the porous medium: ∂p K control Ev − φμ ∂t

∂2 p ∂2 p ∂2 p þ þ ∂x2 ∂y2 ∂z2

! ¼ 0:

ð10Þ

Eq. (10) indicates that the hydraulic diffusivity of fluid pressure in the fluid-saturated porous medium can be expressed in the following form: Dhydraulic ¼

K control Ev : φμ

ð11Þ

If the characteristic length of a fluid flow process in the mineralizing system is defined as Lporefluid − flow, then the characteristic time (i.e. the time-scale) of the fluid flow process can be expressed as follows: 2

t poref luid−f low ¼

2

Lporef luid−f low φμLporef luid−f low ¼ : Dhydraulic K control Ev

ð12Þ

Eq. (12) indicates that as the controlling permeability for different rocks can vary within a range of several orders of magnitude, the characteristic time of a fluid flow process is mainly controlled by the controlling permeability of the porous rocks. To demonstrate this point, it is assumed that the fluid has the following parameters: the bulk modulus of elasticity is 2 × 10 9Pa; the dynamic viscosity is 1 × 10 − 3(Pa ⋅ s) and the porosity of the porous rock is 0.1. For an upper crust of 10 km thick, the characteristic time of a fluid flow process is 5 × 10 8s for a rock with a controlling permeability of 1 × 10 − 14m 2. The characteristic time is reduced to 5 × 10 4s when the rock has a controlling permeability of 1 × 10 − 10m 2. This means that fluid pressure in the pore space of a porous rock diffuses much slower than stress does in the solid skeleton of the porous rock unless the porous rock is much more permeable, such as a well-sorted sandstone. The characteristic time (i.e. the time-scale) of a heat transfer process associated with a mineralizing system in the field of computational hydrodynamics of mineralization can be estimated by the thermal diffusivity of the fluid-saturated porous medium. Through deriving the energy balance equation of a mineralizing system, the thermal diffusivity of a fluid-saturated porous medium due to heat conduction can be expressed as follows: Dconduction ¼ 

λe ρcp



ð13Þ

e

where Dconduction is the thermal diffusion of the fluid-saturated porous medium due to heat conduction; (ρcp)e is a comprehensive parameter to represent the product of the density and specific heat in the fluid-

26

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

saturated porous medium and λe is the equivalent thermal conductivity of the fluid-saturated porous medium. If the thermal equilibrium is reached between the fluid and solid skeleton in the fluidsaturated porous medium (Nield and Bejan, 1992), both (ρcp)e and λe can be expressed as porosity-weighted combinations of the fluid and solid properties: λe ¼ φλf þ ð1−φÞλs ;

  ρcp ¼ φρf cpf þ ð1−φÞρs cps e

ð14Þ

where λf is the thermal conductivity of the fluid in the porous medium; λs is the thermal conductivity of the solid skeleton in the porous medium and φ is the porosity of the porous medium; ρf and ρs are the densities of the fluid and solid skeleton; cpf and cps are the specific heats of the fluid and solid skeleton, respectively. If the characteristic length of a heat transfer process in the mineralizing system is defined as Lthermal, then the characteristic time (i.e. the time-scale) of the heat transfer process due to heat conduction can be expressed as follows:

t heat−conduction

  ρcp L2thermal L2thermal e ¼ ¼ : Dconduction λe

ð15Þ

As a rough approximation for an upper crust of 10 km thick, the following parameters are used to estimate the characteristic time of the heat transfer process due to heat conduction in a fluid-saturated porous medium: the densities of the fluid and solid skeleton are 1000 kg/m 3 and 2600 kg/m 3; the specific heats of the fluid and solid skeleton are 4184J/(kg ⋅ °C) and 878J/(kg ⋅ °C) respectively; the thermal conductivity of the fluid is 0.59 W/(m ⋅ °C); the thermal conductivity of the solid skeleton in the porous medium is 2.9 W/(m ⋅ °C) and the porosity of the porous medium is 0.1. These parameters yield an equivalent thermal conductivity of 2.67W/(m ⋅ °C) for the fluidsaturated porous medium and a characteristic time of 2.475 × 10 14 s, which is equivalent to about 7.79 million years, for the heat transfer process due to heat conduction in the fluid-saturated porous medium. Unlike heat conduction, the characteristic time of the heat transfer process due to heat advection/convection in a mineralizing system is mainly controlled by the average velocity of advective/convective fluid flow in the mineralizing system. This leads to the following equation for estimating the characteristic time of a heat transfer process due to heat advection/convection in a mineralizing system:

t heat−advection=convection

  ρcp Lthermal e  ¼ ρf cpf V advection=convection

ð16Þ

where theat − advection/convection is the characteristic time of a heat transfer process due to heat advection/convection in the mineralizing system; V advection=convection is the average Darcy velocity of either advective or convective fluid flow in the mineralizing system. If the average velocity of convective fluid flow is equal to 5 cm per year (i.e. about 1.58 × 10 − 10m/s), which is of the same order of magnitude as the convective fluid flow resulting from a typical temperature gradient in the crust (Zhao et al., 2008a), then the characteristic time of the heat transfer process due to heat convection is about 3.11 × 10 13s, which is equivalent to about 0.99 million year, for an upper crust of 10 km thick. This demonstrates that in terms of heat transfer, convective/advective heat transfer is much more efficient than conductive heat transfer in a mineralizing system within the upper crust of the Earth. Due to the similarity between a heat transfer equation and a mass transport equation, the characteristic time (i.e. the time-scale) of a mass transport process associated with a mineralizing system in the field of computational hydrodynamics of mineralization can be

estimated using the same procedure as that used for estimating the characteristic time of a heat transfer process in the fluid-saturated porous medium. This means that the solute diffusivity of the porous medium can be used to derive the characteristic time of a mass transport process due to mass (i.e. solute) diffusion in the mineralizing system as follows: 2

t diffusion ¼

Lmass−transport Ddiffusion

ð17Þ

where tdiffusion is the characteristic time of a mass transport process due to mass diffusion in the mineralizing system; Ddiffusion is the solute diffusivity in the porous medium and Lmass − transport is the characteristic length of the mass transport process in the fluid-saturated porous medium. Eq. (17) indicates that if the solute diffusivity in a porous medium is 1 × 10 − 9m 2/s, the characteristic time of a mass transport process due to mass diffusion in a mineral system is about 1 × 10 17s, which is equivalent to 3150 million years, when the characteristic length of the mineral system is 10 km. In the case of considering solute dispersion in the fluid-saturated porous medium, the characteristic time of a mass transport process due to mass dispersion in the mineralizing system can be estimated using the following formula: t dispersion ¼

φL2mass−transport  α V advection=convection

ð18Þ

where tdispersion is the characteristic time of a mass transport process due to mass dispersion in the mineralizing system; α is the solute dispersivity in the porous medium; φ is the porosity of the porous medium; V advection=convection is the average Darcy velocity of either advective or convective fluid flow in the mineralizing system; and Lmass − transport is the characteristic length of the mass transport process in the fluidsaturated porous medium. If the solute dispersivity in the porous medium is equal to 10 m and the average Darcy velocity of convective fluid flow is equal to 5 cm per year (i.e. about 1.58 × 10 − 10m/s), then the characteristic time of a mass transport process due to mass dispersion in a mineralizing system is about 6.33 × 10 15s, which is equivalent to about 199.4 million years, when the characteristic length of the mineral system is 10 km and the porosity of the porous medium is 0.1. Similarly, the characteristic time of a mass transport process due to advective/convective mass transport in a mineralizing system can be expressed as follows: φLmass−transport t mass−advection=convection ¼  V advection=convection

ð19Þ

where tmasst − advection/convection is the characteristic time of a mass transport process due to advective/convective mass transport in the mineralizing system; V advection=convection is the average Darcy velocity of either advective or convective fluid flow in the mineralizing system; and φ is the porosity of the porous medium. To compare the characteristic time of a mass transport process due to convective mass transport with that due to mass dispersion, the same parameters as those used to estimate the characteristic time of a mass transport process due to mass dispersion are employed here. Thus, it is assumed that the average Darcy velocity of convective fluid flow is equal to 1.58 × 10 − 10m/s and that the porosity of the porous medium is 0.1. As a result, the characteristic time of the mass transport process due to convective mass transport in the mineralizing system is about 6.33 × 10 13s, which is equivalent to about 2 million years, when the characteristic length of the mineral system remains 10 km. Since the characteristic time of a mass transport process due to convective mass transport is much smaller than that due

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

to mass dispersion, it has demonstrated that for a mass transport process in a mineralizing system, convective/advective mass transport is much more efficient than that due to mass dispersion to transport dissolved species within the mineralizing system. In terms of the characteristic time (i.e. the time-scale) of a chemical reaction process evolved in a mineralizing system, it can be determined by the chemical reaction rate constant, from the chemical kinetics point of view. Due to the highly variable nature of chemical reaction rates, the characteristic time of a chemical reaction process in a mineralizing system can vary from a few nano-seconds (in the case of a homogeneous reaction) to a few million years (in the case of a heterogeneous reaction). Nevertheless, since the reaction rate of a chemical reaction mainly affects the source/sink terms of reactive mass transport equations, this large range of time scales has little effect on the numerical stability when reactive mass transport equations are solved numerically. 4. Simulation of fluid convection associated with a mineral forming system Convective fluid flow is important for a mineralizing system in the field of computational hydrodynamics of mineralization, at least from the following three points of view. (1) Since the fluid flows in a closed loop within a top-sealed mineralizing system, consumption of the fluid within the mineralizing system is minimal. This enables convective fluid flow to last for relatively long periods of time, provided that elevated heat flow can be maintained at the bottom of the mineralizing system. (2) Since the fluid flows in a closed loop, it is an effective and efficient process for mixing different dissolved species within the mineralizing system. (3) Convective fluid flow may result in highlylocalized temperature anomalies in the mineralizing system. This provides a favorable condition under which high-grade, giant ore deposits may form in the upper crust of the Earth. In terms of investigating convective fluid flow in fluid-saturated porous rocks, much has been done over the past years (Alavyoon, 1993; Bau and Torrance, 1982; Bjorlykke et al., 1988; Caltagirone and Bories, 1985; Chevalier et al., 1999; Gasser and Kazimi, 1976; Horne and Caltagirone, 1980; Horton and Rogers, 1945; Kaviany, 1984; Lapwood, 1948; Lebon and Cloot, 1986; Pillatsis et al., 1987; Tournier et al., 2000; Yang, 2006; Yang et al., 2010; Zhao et al., 1997, 2004, 2006). Nield and Bejan have summarized the early work carried out before 1992 on convective fluid flow in porous media in their book (Nield and Bejan, 1992). Under the stimulus of understanding the physical and chemical processes that control the formation of ore deposits in the upper crust of the Earth, Zhao and his coworkers have carried out extensive and systematic studies on convective and advective heat transfer in geological systems over the past 14 years. As a result, Zhao et al. (2008a) have published a monograph to summarize their research achievements on the theoretical analysis and numerical simulation of fluid flow in geological systems of crustal scales. This large amount of research has been carried out, both theoretically and numerically, to investigate the detailed physical mechanisms related to the condition to trigger the convective fluid flow, and the effects of geological heterogeneities on the convective heat transfer in the crust of the Earth. Generally, the following main conclusions have been made in the monograph of Zhao et al. (2008a): (1) if the Rayleigh number of the system is either critical or supercritical, convective fluid flow can take place in the fluid-saturated porous rocks where the fluid pressure gradient is close to the hydrostatic pressure gradient. (2) If the fluid pressure gradient is close to the lithostatic pressure gradient, convective fluid flow cannot take place in the porous medium that has constant (but different values of) temperature and impermeable boundary conditions at both the top and the bottom of the system, but it can take place if the system has a permeable top with the constant pressure and temperature, and a

27

bottom with the constant upward fluid velocity and conductive heat flux. (3) The heterogeneity of either the medium permeability or the thermal conductivity has a significant effect on the convective heat transfer within the system. The material thermoelasticity may also affect the heat transfer, depending on the relative hardness of the rock masses. (4) Convective fluid flow can also take place in the geological fault zone, depending on the critical Rayleigh number of the system, which is a function of the ratio of the fault height to the fault thickness. From the physical point of view, the occurrence of convective fluid flow in the Earth's crust is an emerging behavior of the complex Earth system. Once the Rayleigh number reaches the critical Rayleigh number of the system, convective fluid flow takes place instantaneously. The key point associated with this emerging phenomenon is that the resulting convective fluid flow pattern consistent with the critical Rayleigh number is of a steady-state nature and therefore, is independent of time if the thermodynamic properties of the system remain unchanged. This indicates that convective fluid flow can be analyzed mathematically and simulated numerically in a steady-state analysis of the system. Based on this recognition, the steady-state approach was used to theoretically investigate the convective instability of fluid flow both in a block model of the Earth's crust and in three-dimensional geological fault zones (Zhao et al., 2008a). To use the steady-state approach efficiently, Zhao et al. (1997) developed the progressive asymptotic approach procedure to numerically simulate the steady-state convective fluid flow in the Earth's crust which is represented as fluid-saturated porous media. Since convective fluid flow is an intrinsic feature of the complex Earth system, it can be found through different kinds of perturbation methods. In this regard, one can use the transient-state approach to numerically simulate the system as a transient system with an arbitrary perturbation of temperature. If the numerical algorithm is robust enough, convective fluid flow should take place once the Rayleigh number is equal to or greater than the corresponding critical Rayleigh number of the system. In this case, the time-varying fluid flow process is an artifact associated with the transient-state approach and therefore, has no true physical meaning. For example, one can use different temperature perturbations in the transient numerical simulation, so that one could obtain different time-varying fluid flow processes, even if the system approaches the same convective fluid flow pattern as predicted theoretically. Therefore, the real value of the transientstate approach to the convective instability of fluid flow is to find the critical Rayleigh number and the related convective fluid flow pattern. Nevertheless, compared with the steady-state approach, the transient-state approach is computationally inefficient as a result of using many time steps in the numerical simulation. Thus, if one is only interested in obtaining the convective fluid flow pattern, then the steady-state approach is superior to the transient-state approach for the computational simulation of convective fluid flow in the upper crust of the Earth. Using the steady-state approach and the progressive asymptotic approach algorithm, Zhao et al. (2008a) simulated the convective fluid flow within an upper crust of 100 km long and 13 km thick. The computational domain is modeled by 3956 quadrilateral elements with the following boundary conditions: (1) temperature at the top and the bottom of the computational domain is 20 °C and 150 °C respectively; (2) all the lateral vertical boundaries are insulated and impermeable in the horizontal direction; and (3) both the top and bottom boundaries are impermeable in the vertical direction. The following parameters are used in the computational simulation: the dynamic viscosity of the fluid is 10 − 3N ⋅ s/m 2; the volumetric thermal expansion coefficient of the fluid is 2.07 × 10 − 4(1/°C); the specific heat of the fluid is 4185 J/(kg ⋅ °C); the thermal conductivity coefficients of the fluid in both the horizontal and vertical directions are 0.6 W/(m ⋅ °C); the porosity of

28

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

(Fluid velocity)

(Temperature) Fig. 1. Distributions of fluid flow and temperature in the crustal model: in the case of showing fluid velocity, red color represents high velocity regions, while blue represents low velocity regions; in the case of showing temperature, red color represents high temperature regions, while blue represents low temperature regions.

the porous rock matrix is 0.1; the specific heat of the porous rock matrix is 815 J/(kg ⋅ °C); the thermal conductivity coefficients of the porous rock matrix in both the horizontal and vertical directions are 3.35 W/(m ⋅ °C); the permeability of the porous rock is 10 − 14m 2. Fig. 1 shows the distributions of the fluid velocity and temperature in the crustal model, where the maximum fluid velocity is 3.7 × 10 − 10m/s and the maximum temperature is 150 °C. In this figure, the unit of fluid velocity is m/s, while the unit of temperature is °C. From the numerical results, it can be observed that: (1) multiple convection cells can be reasonably simulated in such a crustal model, and (2) the distribution of temperature is highly localized. This is a direct result of the coupling between the temperature and the fluid velocity within the convective fluid flow. From the modern mineralization theory (Zhao et al., 2008a), the rate of fluid–rock mass transfer is directly proportional to the scalar product of the fluid flow velocity and temperature gradient. This means that the peak region (where the fluid flow is upward) of the temperature distribution is a favorable location for mineralization in the computational model if the variation of mineral solubility with temperature (and pressure) under consideration is prograde (i.e. the mineral solubility increases with the increase of temperature (and pressure)). Due to the important roles that are played by convective fluid flow in mineralizing systems, a considerable amount of research has been carried out to describe mineralization patterns in mineralizing systems involving convective fluid flow regimes (Fu et al., 2010; Harcouet-Menou et al., 2009; Magri et al., 2010; Sorjonen-Ward et al., 2002; Yang et al., 2006; Zhao et al., 1998, 2008b). In particular, the related computational simulation results have demonstrated that convective fluid flow is the dynamic mechanism of resulting equal-distance distribution of Au deposits within the Bardoc fault zone of the Yilgarn Craton, Western Australia (Zhao et al., 2008b). Although more complicated reactive mass transport processes are considered in three-dimensional mineralizing systems involving convective fluid flow regimes (Alt-Epping and Zhao, 2010), the feedback effect of mineral precipitation/dissolution on the convective fluid flow pattern was neglected in previous computational models. Since mineral precipitation/dissolution causes a change in the porosity of a porous medium, this change can result in a change in the microscopic flow channels, so that the convective fluid flow pattern may change as the mineral precipitation/dissolution proceeds in the mineralizing system. Hence, in order to accurately predict mineralization patterns in convective fluid flow systems, it is necessary to consider the feedback effect of mineral precipitation/dissolution on the convective fluid flow pattern in the mineralizing system of a convective fluid flow regime.

5. Simulation of fluid focusing and mixing associated with a mineral forming system Structure controls on mineralization in the upper crust of the Earth have been an important topic for many years (Cox, 1999; Hodgson, 1989; Leach and Rowan, 1986; Matthai, 2003; Norton and Knight, 1977; Oliver, 1995; Phillips, 1972; Valenta et al., 1994; Zhang et al., 2010). Since the main focus of this topic is on the fluid pressure gradient driven flow, the mechanical deformation of crustal rocks needs to be considered. Generally, the fluid pressure gradient driven flow is closely associated with overpressure in the crust and can be caused by the following mechanism: (1) tectonic deformation (Leach and Rowan, 1986; Oliver, 1986, 1995; Sibson, 1987; Sibson et al., 1975); (2) topographic change (Deming and Nunn, 1991; Garven, 1995) and (3) metamorphic reactions of the rock (Cathles, 1981; Cathles and Smith, 1983; Ge and Garven, 1992). One of the key factors affecting the fluid pressure gradient driven flow is the distribution and generation of flow channels in the crustal rocks. Since permeable fault and/or fracture zones have enormous capability to create efficient fluid pathways in the crust, they are often the most important structures to control mineralization in many ore deposits (Harcouet-Menou et al., 2009; Matthai, 2003; Sibson, 1994; Sibson et al., 1988; Yang, 2006; Yang et al., 2010; Zhang et al., 2003). Except for fluid flow focusing in permeable porous rocks, the mixing of two or more fluids is commonly suggested as a mechanism for precipitating minerals from solution in the permeable porous rocks. Examples include Uranium deposits (Wilde and Wall, 1987), MVT deposits (Appold and Garven, 2000; Bethke, 1986; Garven et al., 1999), Irish Pb-Zn deposits (Everett et al., 1999; Hitzman, 1995), Vein Gold deposits (Cox et al., 1995; Matthai et al., 1995) and Carlin Gold deposits (Cline and Hofstra, 2000). The mixing process is attractive because it enables two fluids of contrasting fluid compositions to mix and hence generate chemical conditions favorable for mineral precipitation. For a given fault (either permeable or impermeable relative to its surrounding rock) of elliptic shape in a hydrodynamic system, a parameter, known as the flow focusing factor, is used to measure the flow contrast between the fault and its surrounding rock. The flow focusing factor is defined as the ratio of the maximum fluid flow velocity in the fault to the inflow velocity within the surrounding rock (Zhao et al., 2008a). Fig. 2 shows the variation of the flow focusing factor with the aspect ratio due to different permeability ratios when the inflow is parallel to the long axis of the elliptic fault. In this figure, the permeability ratio (represented by PR) is defined as the ratio of the fault permeability to the surrounding rock permeability, while the aspect ratio is

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

Log(Flow Focusing Factor)

4

faults, even if the load applied to the system is large enough to cause fault propagation. To overcome this weakness, it is necessary to develop advanced numerical methods and algorithms to effectively and efficiently simulate random crack generation and propagation in brittle rocks (Zhao et al., 2008c, 2009).

PR=5 PR=10 PR=100 PR=1000 PR=10000

3

2

6. Simulation of reactive mass transport associated with a mineral forming system

1

Due to broad involvement in many scientific and engineering fields, the research work associated with the reactive mass transport modeling has been extensively carried out over the past years (AltEpping and Zhao, 2010; Appold and Garven, 2000; Garven and Freeze, 1984; Lichtner, 1996; Matthai et al., 2010; Mills et al., 2007, 2009; Raffensperger and Garven, 1995; Steefel and Lasaga, 1994; Steefel and Yabusaki, 1996; Yeh and Tripathi, 1991; Zhao et al., 1998, 2001, 2006, 2007a). Certainly, the review of this subject alone deserves another paper. The focus of this section is selectively on the following three emerging issues related to the reactive mass transport modeling of mineralizing systems: (1) modeling of the chemical dissolution-front instability in a mineral dissolution system; (2) implementation of new equations of state and constitutive relationships as well as the dependence of thermodynamic parameters, such as density, viscosity, diffusivity, and thermal conductivity, on the temperature, pressure and chemical environment of a mineralizing system; and (3) simultaneous modeling of the multi-scales and multi-processes problems associated with a mineralizing system. Most mineralizing systems in the field of computational hydrodynamics of mineralization may involve mineral dissolution/precipitation reactions. In the case of a mineral dissolution reaction, the chemical dissolution front in a fluid-saturated porous medium can become unstable if the geochemical system under consideration is supercritical, in which a simple planar dissolution front can involve into fingering shapes, resulting in flow focusing in the chemical dissolution system. Since flow focusing plays an important role in the formation and enrichment of large ore deposits, the instability of a chemical dissolution front may become a potential mechanism for ore body formation and mineralization in the upper crust of the Earth. This scientific problem is often called the chemical dissolution-front instability problem and has been studied for several years (Zhao et al., 2008d, 2009). To represent the fundamental characteristics of a chemical-dissolution reaction system, the following Zhao number was proposed:

0 1

2

3

4

5

Log(Aspect Ratio) Fig. 2. Variation of flow focusing factor with aspect ratio due to different permeability ratios (PR): inflow parallel to the long axis of the inclusion.

defined as the ratio of the fault length to its thickness. From this result, it can be concluded that for a given aspect ratio of the fault, the flow focusing factor increases with an increase in the permeability ratio of the fault to the surrounding rock until it reaches the corresponding limiting value. Similarly, for a given permeability ratio of the fault to the surrounding rock, the flow focusing factor increases with an increase in the aspect ratio of the fault until it reaches its corresponding limiting value. To investigate the potential mineralization patterns within a permeable fault and its surrounding rock, Zhao et al. (2007a, 2009) proposed a generic fault model, in which two fluids were injected from the two sides at one tip of the fault. In this generic model, the interaction between solute advection, solute dispersion/diffusion and chemical kinetics was considered in an analytical manner. Based on three time scales of the three processes involved in the model, the following dimensionless parameters are defined: φkR l V ¼ Time Scale f or Solute Advection=Time Scale f or Chemical Reaction

Da ¼

ð20Þ Z¼

29

k R l2 ¼ ðTime Scale f or Solute Dispersion=Dif f usionÞ= D Time Scale f or Chemical Reaction

ð21Þ

where Da is the Damköhler number; V is the characteristic Darcy velocity of the system; D is the solute diffusion/dispersion coefficient; l is the characteristic length of the controlling process in the system; kR is the reaction rate with a unit of [s − 1]; and ϕ is the porosity of the porous medium. Using these two dimensionless parameters, it is possible to investigate the potential mineralization patterns within the permeable fault and its surrounding rock. As shown in Fig. 3, there are three possible cases in which mineral saturation can be attained. As a result, three types of mineralization patterns can be described as follows: (1) in the case of type 1 mineralization pattern, mineral precipitation comprises a thin lenticular shape within the fault zone and starting at the lower tip of the fault; (2) in the case of type 2 mineralization pattern, mineral precipitation cannot take place within the permeable fault zone; and (3) in the case of type 3, mineral precipitation has a considerable thickness within the fault zone. Because of a significant size difference between the fault length and thickness, the adaptive mesh refinement scheme needs to be used to simulate fluid flow, heat transfer and mass transport within a fault. Since current modeling practice is to simulate faults with described geometric shapes, it does not consider the propagation of the

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u V p vflow t ffi Zh ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   k  chemical A p C eq φf D φf

ð22Þ

where vflow is the Darcy velocity of the fluid flow; φf and D(φf) are the final porosity of the porous medium and the corresponding diffusivity of the chemical species after the completion of dissolvable-mineral dissolution; Ceq is the equilibrium concentration of the chemical species; V p is the average volume of the dissolvable minerals; Āp is the averaged surface area of the dissolvable minerals; and kchemical is the rate constant of the controlling chemical dissolution reaction. For the purpose of understanding the physical meanings of each term in the Zhao number, Eq. (22) can be rewritten in the following form: Zh ¼ F Advection F Diffusion F Chemical F Shape

ð23Þ

where FAdvection is a term to represent the solute advection; FDiffusion is a term to represent the solute diffusion/dispersion; FChemical is a term to represent the kinetics of the chemical dissolution reaction; and FShape is a term to represent the shape factor of the dissolvable

30

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

Da = Time Scale for Solute Advection / Time Scale for Chemical Reaction

Solute Diffusion Dominates

Chemical Reaction Dominates

Type 1

Type 1

Chemical Reaction Dominates

Type 3

Type 2

1.0 Type 2

Solute Advection Dominates

Type 2

Type 3 0 1.0

Z = Time Scale for Solute Diffusion / Time Scale for Chemical Reaction Fig. 3. Three types of mineralization patterns in a vertical fault zone due to the interaction between solution advection, solution dispersion/diffusion and chemical kinetics: red color represents mineral precipitation regions, and dashed lines represent the locations of the fault zone.

mineral in the chemical dissolution reaction system. These terms can be expressed as follows: F Advection ¼ vflow

ð24Þ

1 F Diffusion ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi φ f D φf

ð25Þ

F Chemical

F Shape

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ¼ kchemical C eq

ð26Þ

vffiffiffiffiffiffi u uV p ¼ t : A

ð27Þ

p

Eqs. (24) to (27) indicate that the Zhao number is a dimensionless number that can be used to represent the geometrical, hydrodynamic, thermodynamic and chemical kinetic characteristics of a chemical dissolution reaction system in a comprehensive manner. Note that the critical condition, under which a chemical dissolution front becomes unstable, can be expressed by the critical value of this dimensionless number.

Zhcritical ¼

ð3−βÞð1 þ βÞ ; 2ð1−βÞ

kðφ Þ β ¼  0 k φf

ð28Þ

where Zhcritical is the critical Zhao number of the chemical dissolution reaction system; k(φ0) is the initial permeability corresponding to the initial porosity of the porous rock; k(φf) is the final permeability corresponding to the final porosity, φf, of the porous rock. Both the Zhao number and the corresponding critical one can be used to assess the instability likelihood of a chemical dissolution front involved in the mineralizing system within the upper crust of the Earth. If the Zhao number of a mineralizing system is smaller than its corresponding critical Zhao number, then the mineralizing system is subcritical. Otherwise, it is either critical or supercritical. Note that although the stability criterion, which is used to assess the instability likelihood of the chemical dissolution front in a geochemical system, is established from considering a single chemical dissolution reaction in the geochemical systems, it can be approximately extended to the assessment of chemical dissolution-front instability in a geochemical system involving several mineral dissolution reactions. From

this point of view, a controlling mineral dissolution reaction in a mineralizing system needs to be identified so that the established stability criterion can be used to judge the instability of the chemical dissolution front associated with this particular controlling mineral dissolution reaction in the mineralizing system. On the other hand, the previous studies (Zhao et al., 2009, 2010) have demonstrated that the porosity enhancement patterns obtained from a subcritical dissolution system are significantly different from those obtained from a supercritical dissolution system. Since conventional numerical algorithms are usually developed for simulating reactive mass transport problems involving subcritical chemical dissolution reactions in a geochemical system (Appold and Garven, 2000; Garven and Freeze, 1984; Hammond and Lichtner, 2010; Hammond et al., 2011; Lichtner, 1996; Matthai et al., 2010; Mills et al., 2007, 2009; Raffensperger and Garven, 1995; Steefel and Lasaga, 1994; Steefel and Yabusaki, 1996; Yeh and Tripathi, 1991), they are well capable of simulating complex dissolution/precipitation chemical reactions and the resulting porosity/permeability changes for a certain range of temperature, pressure and ionic strengths in a subcritical chemical dissolution system, but may be unsuitable to simulate the instability behavior of the chemical dissolution front in a chemical dissolution dominated system of the supercritical Zhao number. Thus, some special numerical algorithms may be needed to appropriately simulate the morphological evolutions of chemical dissolution fronts in such supercritical dissolution systems. The phenomenon of dissolution front instability has been not only observed in simple geochemical systems, such as the infiltration of a low-pH fluid into a carbonate rock (e.g. the Karst formation), but also successfully applied to create porosity and preferential flow channels to increase oil production in petroleum industry through injecting acid into carbonate rocks (Cohen et al., 2008; Kalia and Balakotaiahand, 2007). However, since ore forming systems may involve a set of multiple parallel dissolution/precipitation reactions, which can result in a porosity/permeability increase due to mineral dissolution as well as a porosity/permeability decrease due to mineral precipitation, the phenomenon of dissolution front instability might be swallowed in some ore forming systems unless dissolution reactions are predominant processes during a certain period of the ore forming system. When hydrothermal mineralizing systems in the upper crust of 10 km are considered, the maximum temperature and (hydrostatic) pressure may reach about 300°C and 1 kba, respectively. However, when magmatic mineralizing systems involving the whole crust of 50 km are taken into account, the maximum temperature and (hydrostatic) pressure may reach about 1500°C and 5 kba, respectively. This means that the variations of fluid properties, which are commonly temperature and pressure dependent, may have significant effects on the thermodynamic and chemical behaviors of a mineralizing

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

system. Although numerical models have been widely used to simulate hydrothermal mineralizing systems in the upper crust, the numerical model (even of generic nature) that can be used to simulate magmatic mineralizing systems in the whole crust remains very limited. Since any numerical model involves fluid properties and thermodynamic data, another important fundamental issue that researchers involved in modeling hydrothermal or mineralized systems are presently concerned with is to realistically represent fluid properties and chemical reactions at high pressures, high temperatures and high ionic strengths as one would find in ore forming systems. While numerous efforts have been made to develop new equations of state and new fluid property models in recent years (Barnes, 1997; Geiger et al., 2006a,b; Pokrovskii, 1999; Shock and Helgeson, 1988; Shock and Koretsky, 1995; Wang et al., 2010), some are with focuses on systems of geochemical importance such as aqueous solutions containing ions or neutral molecules that exist in significant amount in mineralizing systems (Anderko and Lencka, 1998; Lencka et al., 1997; Wang and Anderko, 2008; Wang et al., 2002). Such models can provide not only various fluid properties that are relevant to modeling hydrodynamics of mineralization, such as density, viscosity, surface tension, diffusivity, and thermal conductivity, but also the thermodynamic properties of chemical reactions that are necessary for describing the fluid speciation and mineral dissolution/precipitation reactions, as well as the energy/heat absorbed or released from chemical reactions. To simulate a mineralizing system more realistically, such new equations of state and new fluid property models need to be implemented in the computational simulation of the mineralizing system so that fluid properties and chemical reactions at high pressures, high temperatures and high ionic strengths can be realistically simulated in ore forming systems. Because of significant differences between engineering problems and mineralization problems, commercial computer codes and related algorithms, primarily designed for solving engineering problems, cannot be directly used to solve mineralization problems without some modifications. For this reason, it is necessary either to develop new computer codes for solving multi-scales and multi-processes mineralization problems or to modify the existing commercial computer codes, originally designed for solving engineering problems, so as to be suitable for solving mineralization problems. For example, Steefel and Yabusaki (1996) developed a computer code known as OS3D to solve the coupled problems between fluid flow, heat transfer and reactive transport processes that occur in ore forming systems through coupling flow solvers with external geochemistry solvers. Matthai et al. (2004) developed a computer code known as the complex systems platform (CSP) to accurately and efficiently model multiphase flow in mineralizing systems, on the basis of a fully conservative node-centered finite volume method coupled with the Galerkin finite element method on an unstructured triangular grid with a complementary finite volume subgrid. They also employed the non-uniform rational B-splines (NURBS), which is an advanced technique in the field of computing geometry, to simulate complicated geometric shapes in ore forming systems and other geological systems (Paluszny et al., 2007). In recent years, Mills et al. (2007, 2009) and Hammond and Lichtner (2010, 2011) developed a computer code known as PFLOTRAN for modeling multi-phase, multi-component flow and reactive transport through using massively parallel computers. PFLOTRAN is built on top of the Portable, Extensible Toolkit for Scientific Computation (PETSc) parallel platform so that it exhibits excellent performance on the largest-scale supercomputers. In addition, PFLOTRAN also employs domain-decomposition parallelism, which is a modern technique available for large-scale scientific and engineering computation in the field of computational science. This means that even though PFLOTRAN is still under development and has not been applied to the simulation of ore-forming systems, it is computationally capable of simulating large-scale ore forming systems in the future. Although many research computer codes have

31

been developed to deal with some aspects of multi-scale and multiprocess mineralization problems (Bethke et al., 1993; Diersch, 2002; Hammond and Lichtner, 2010; Hammond et al., 2011; Matthai et al., 2004; Mills et al., 2007, 2009; Molson and Frind, 1993; Pruess, 1992; Steefel and Yabusaki, 1996), there is not any one code available that can be used to simultaneously simulate the fully coupled problem between material deformation, fluid flow, heat transfer, mass transport and chemical reactions in a mineralizing system involving all multi-scales and multi-processes. From the computational simulation point of view, the current status of reactive mass transport simulation is based on the following two main approaches. In the first approach, mass transport equations (i.e. partial differential equations), in which solute advection and diffusion/dispersion are considered but the source/sink terms due to chemical reactions are neglected, are firstly solved. To consider the effects of the source/sink terms due to chemical reactions to the reactive mass transport process, chemical reaction equations are then solved. As a result, the mass transport and the chemical reaction are considered separately and sequentially in the computational simulation of the reactive mass transport process (Alt-Epping and Smith, 2001; Steefel and Yabusaki, 1996). In the first approach, the originally coupled reactive-mass-transport processes are artificially divided into two separate steps. In the second approach, the source/sink terms due to chemical reactions are explicitly included in mass transport equations, so that both mass transport and chemical reactions are considered simultaneously in the computational simulation of the reactive mass transport process (Zhao et al., 2001). The major advantage of using the first approach is that fluid flow solvers can be easily coupled with external geochemistry solvers (e.g. EQ3/EQ6 (Wolery, 1979)) through a common interface. The major disadvantage is that nonlinear interaction between mass transport and chemical reactions can only be approximated. This approximation implies that mass conservation cannot be strictly satisfied in the computational simulation of the reactive mass transport process. In contrast, the second approach allows the direct consideration of nonlinear interaction between mass transport and chemical reactions in a seamless manner (Zhao et al., 2009). The major disadvantage of using the second approach is that a source/sink term must be explicitly included in each of the reactive mass transport equations. Nevertheless, in order to ensure the correctness of computational simulation results for reactive mass transport processes, it is necessary to further develop software packages with more functionality to solve reactive mass transport equations simultaneously. This should become a priority research direction in the field of computational hydrodynamics of mineralization and include the following aspects: (1) developing advanced/robust simultaneous solvers and adaptive mesh refinement schemes so as to create a computer platform that can be used to simultaneously solve the fully coupled problem between material deformation, fluid flow, heat transfer, mass transport and chemical reactions in a mineralizing system; (2) implementing new equations of state and constitutive relationships as well as the dependence of thermodynamic parameters on the temperature, pressure and chemical environment of a mineralizing system so as to simulate a mineralizing system more realistically; (3) developing advanced numerical methods and algorithms to simulate random crack generation and propagation in brittle rocks so that magma intrusion processes associated with mineralizing systems can be appropriately simulated; and (4) developing advanced mesh generators to make it possible to efficiently simulate the complicated three-dimensional geometry and complex material distribution of a mineralizing system within the upper crust of the Earth. 7. Conclusions and discussions Computational hydrodynamics of mineralization is an important sub-discipline of the emerging computational geoscience. Under

32

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

the stimulus of exploring concealed mineral resources in the deep Earth, this sub-discipline has received much attention and considerable progress has been made in the recent years. The main purpose of establishing computational hydrodynamics of mineralization is to simulate the contributions of fluid-flow related processes to the formation of ore deposits within the upper crust of the Earth. The key aspects of research in the field of computational hydrodynamics of mineralization is to understand, both qualitatively and quantitatively, how, why and where fluid flow can play important roles in the formation of concealed giant mineral deposits in the deep Earth. This may establish some theoretical basis for finding giant ore deposits through geochemical exploration and other methods. Due to the complex and complicated nature of a mineralizing system, the following fundamental issues should be carefully considered in the computational simulation of the mineralizing system. First, multiple roles played by fluids indicate that ore body formation in a mineralizing system is the direct consequence of diverse processes such as material deformation, fluid flow, heat transfer, mass transport and chemical reactions. Second, as each of these processes can involve different time and length scales, the mineral system considered in the field of computational hydrodynamics of mineralization is characterized by multiple time and length scales. Third, convective fluid flow is an effective and efficient mechanism to enhance ore body formation and mineralization as a result of the fluid focusing and mixing in the mineralizing system. However, due to mineral dissolution/precipitation, the flow channels of convective fluid are a direct consequence of pore space reorganization in the fluid-saturated porous medium. Fourth, in order to transport minerals from the deep Earth to the shallow (upper) surface of the Earth, reactive mass transport must be appropriately simulated in the related computational models. In this particular aspect, soluble minerals can be added to or removed from the fluid through precipitation/dissolution chemical reactions. Although significant progress has been made on the numerical modeling of mineralizing systems, research on some aspects, such as simultaneous simulation of multi-processes and multi-scales of mineralizing systems and accurate modeling of reactive mass transport processes in mineralizing systems, should be strengthened in the future work. On the other hand, for the purpose of ensuring the correctness of computational simulation results for reactive mass transport processes, it is imperative to further develop software packages with more functionality to simultaneously solve reactive mass transport equations. This should become a priority research direction in the field of computational hydrodynamics of mineralization and include the following aspects: (1) developing advanced/robust simultaneous solvers and adaptive mesh refinement schemes so as to create a computer platform that can be used to simultaneously solve the fully coupled problem between material deformation, fluid flow, heat transfer, mass transport and chemical reactions in a mineralizing system; (2) implementing new equations of state and constitutive relationships as well as the dependence of thermodynamic parameters on the temperature, pressure and chemical environment of a mineralizing system so as to simulate a mineralizing system more realistically; (3) developing advanced numerical methods and algorithms to simulate random crack generation and propagation in brittle rocks so that magma intrusion processes associated with mineralizing systems can be appropriately simulated; and (4) developing advanced mesh generators to make it possible to efficiently simulate the complicated three-dimensional geometry and complex material distribution of a mineralizing system within the upper crust of the Earth. Acknowledgments This work is financially supported by the Natural Science Foundation of China (Grant No: 10872219) and the Western Australian Geothermal

Centre of Excellence. The authors express their thanks to the anonymous referees for their valuable comments, which led to a significant improvement over an early version of the paper.

References Alavyoon, F., 1993. On natural convection in vertical porous enclosures due to prescribed fluxes of heat and mass at the vertical boundaries. International Journal of Heat and Mass Transfer 36, 2479–2498. Alt-Epping, P., Smith, L., 2001. Computing geochemical mass transfer and water/rock ratios in submarine hydrothermal systems: implications for estimating the vigour of convection. Geofluids 1, 163–181. Alt-Epping, P., Zhao, C., 2010. Reactive mass transport modeling of a three-dimensional vertical fault zone with a finger-like convective flow regime. Journal of Geochemical Exploration 106, 8–23. Anderko, A., Lencka, M.M., 1998. Modeling self-diffusion in multicomponent aqueous electrolyte systems in wide concentration ranges. Industrial and Engineering Chemistry Research 37, 2878–2888. Appold, M.S., Garven, G., 2000. Reactive flow models of ore formation in the Southeast Missouri District. Economic Geology and the Bulletin of the Society of Economic Geologists 95, 1605–1626. Barnes, H.L., 1997. Geochemistry of Hydrothermal Ore Deposits. John Wiley and Sons, New York. Bau, H.H., Torrance, K.E., 1982. Low Rayleigh number thermal convection in a vertical cylinder filled with porous materials and heated from below. ASME Journal of Heat Transfer 104, 166–172. Bergins, C., Crone, S., Strauss, K., 2005. Multiphase flow in homogeneous porous media with phase change, Part II: analytical solutions and experimental verification for constant pressure steam injection. Transport in Porous Media 60, 275–300. Bethke, C.M., 1985. A numerical model of compaction-driven groundwater flow and heat transfer and its application to paleohydrology of intracratonic sedimentary basins. Journal of Geophysical Research 90, 6817–6828. Bethke, C.M., 1986. Hydrologic constraints on genesis of the Upper Mississippi Valley Mineral District from Illinois Basin brines. Economic Geology 81, 233–249. Bethke, C.M., Lee, M.K., Quinodoz, H., Kreiling, W.N., 1993. Basin Modeling with Basin2, a Guide to Using Basin2, B2plot, B2video, and B2view. University of Illinois, Urbana. Bjorlykke, K., Mo, A., Palm, E., 1988. Modelling of thermal convection in sedimentary basins and its relevance to diagenetic reactions. Marine and Petroleum Geology 5, 338–351. Bons, P.D., Page, J.D., Elburg, M.A., 2001. Stepwise accumulation and ascent of magmas. Journal of Metamorphic Geology 19, 627–633. Burbidge, D.R., Braun, J., 2002. Numerical models of the evolution of accretionary wedges and fold-and-thrust belts using the distinct-element method. Geophysical Journal International 148, 542–561. Caltagirone, J.P., Bories, S., 1985. Solutions and stability criteria of natural convective flow in an inclined porous layer. Journal of Fluid Mechanics 155, 267–287. Cathles, L.M., 1981. Fluid flow and genesis of hydrothermal ore deposits. Economic Geology, 75th Anniversary Volume, pp. 424–457. Cathles, L.M., Smith, A.T., 1983. Thermal constraints on the formation of Mississippi Valley-type lead-zinc deposits and their implications for episodic basin dewatering and deposit genesis. Economic Geology 78, 983–1002. Chevalier, S., Bernard, D., Joly, N., 1999. Natural convection in a porous layer bounded by impervious domains: from numerical approaches to experimental realization. International Journal of Heat and Mass Transfer 42, 581–597. Clauser, C., 2009. Heat transport processes in the Earth's crust. Surveys in Geophysics 30 (3), 163–191. Cline, J.S., Hofstra, A.A., 2000. Ore-fluid evolution at the Getchell Carline-type gold deposit, Nevada, USA. European Journal of Mineralogy 12, 195–212. Cohen, C.E., Ding, D., Quintard, M., Bazin, B., 2008. From pore scale to wellbore scale: impact of geometry on wormhole growth in carbonate acidization. Chemical Engineering Science 63, 3088–3099. Cox, S.F., 1999. Deformational controls on the dynamics of fluid flow in mesothermal gold systems. Geological Society of London Special Publications 155, 123–139. Cox, S.F., Sun, S.S., Etheridge, M.A., Wall, V.J., Potter, T.F., 1995. Structural and geochemical controls on the development of turbidite-hosted gold quartz vein deposits, Wattle Gully mine, central Victoria, Australia. Economic Geology and the Bulletin of the Society of Economic Geologists 90, 1722–1746. Crone, S., Bergins, C., Strauss, K., 2002. Multiphase flow in homogeneous porous media with phase change, Part I: numerical model. Transport in Porous Media 49, 291–321. Dagan, G., 1989. Flow and Transport in Porous Formations. Springer Verlag, Berlin. Deming, D., Nunn, J.A., 1991. Numerical simulations of brine migration by topographically driven recharge. Journal of Geophysical Research 96, 2485–2499. Diersch, H.J.G., 2002. FEFLOW Reference Manual. Wasy GmbH, Berlin. Donze, F., Mora, P. And, Magnier, S.A., 1994. Numerical simulation of faults and shear zones. Geophysical Journal International 116, 46–52. Driesner, T., Geiger, S., 2007. Numerical simulation of multiphase fluid flow in hydrothermal systems. Reviews in Mineralogy and Geochemistry 65, 187–212. Everett, C.E., Wilkinson, J.J., Rye, D.M., 1999. Fracture-controlled fluid flow in the Lower Palaeozoic basement rocks of Ireland: implications for the genesis of Irish-type Zn–Pb deposits. In: McCaffrey, K.J.W., Lonergan, L., Wilkinson, J.J. (Eds.), Fractures, fluid flow and mineralization: Geological Society, London, Special Publications, 155, pp. 247–276.

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34 Finch, E., Hardy, S., Gawthorpe, R., 2004. Discrete element modeling of extensional fault-propagation folding above rigid basement fault blocks. Basin Research 16, 489–506. Frind, E.O., 1982. Simulation of long-term transient density-dependent transport in groundwater. Advances in Water Resoures 5, 73–88. Fu, F.Q., McInnes, B.I.A., Evans, N.J., Davies, P.J., 2010. Numerical modeling of magmatic– hydrothermal systems constrained by U–Th–Pb–He time–temperature histories. Journal of Geochemical Exploration 106, 90–109. Garven, G., 1995. Continental-scale groundwater flow and geologic processes. Annual Review of Earth and Planetary Sciences 23, 89–117. Garven, G., Appold, M.S., Topygina, V.I., Hazlett, T.J., 1999. Hydrologic modeling of the genesis of carbonate-hosted lead-zinc ores. Hydrogeology Journal 7, 108–126. Garven, G., Freeze, R.A., 1984. Theoretical analysis of the role of groundwater flow in the genesis of stratabound ore deposits: mathematical and numerical model. American Journal of Science 284, 1085–1124. Gasser, R.D., Kazimi, M.S., 1976. Onset of convection in a porous medium with internal heat generation. ASME Journal of Heat Transfer 98, 49–54. Ge, S., Garven, G., 1992. Hydromechanical modeling of tectonically driven groundwater flow with application to the Arkoma foreland basin. Journal of Geophysical Research 97, 9119–9144. Geiger, S., Driesner, T., Heinrich, C.A., Matthai, S.K., 2006a. Multiphase thermohaline convection in the Earth's crust: I. A new finite element–finite volume solution technique combined with a new equation of state for NaCl–H2O. Transport in Porous Media 63, 399–434. Geiger, S., Driesner, T., Heinrich, C.A., Matthai, S.K., 2006b. Multiphase thermohaline convection in the Earth's crust: II. Benchmarking and application of a finite element–finite volume solution technique with a NaCl–H2O equation of state. Transport in Porous Media 63, 417–443. Gerhart, P.M., Gross, R.J., Hochstein, J.I., 1992. Fundamentals of Fluid Mechanics. AddisonWesley Publishing Company, New York. Gow, P., Upton, P., Zhao, C., Hill, K., 2002. Copper-Gold mineralization in the New Guinea: numerical modeling of collision, fluid flow and intrusion-related hydrothermal systems. Australian Journal of Earth Sciences 49, 753–771. Hammond, G.E., Lichtner, P.C., 2010. Field-scale model for the natural attenuation of Uranium at the Hanford 300 area using high performance computing. Water Resources Research 46, 1–31. Hammond, G.E., Lichtner, P.C., Rockhold, M.L., 2011. Stochastic simulation of Uranium migration at the Hanford 300 area. Journal of Contaminant Hydrology 120–121, 115–128. Harcouet-Menou, V., Guillou-Frottier, L., Bonneville, A., Adler, P.M., Mourzenko, V., 2009. Hydrothermal convection in and around mineralized fault zones: insights from two- and three-dimensional numerical modeling applied to the Ashanti belt, Ghana. Geofluids 9, 116–137. Hitzman, M.W., 1995. Mineralisation in the Irish Zn–Pb–(Ba–Ag) Orefield. In: Anderson, K., Ashton, J., Earls, G., Hitzman, M., Tear, S. (Eds.), Irish Carbonate-hosted Deposits, Society of Economic Geologists: Guidebook Series, 21, pp. 25–61. Hobbs, B.E., 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, 505–509. Hodgson, C.J., 1989. The structure of shear-related, vein-type gold deposits: a review. Ore Geology Review 4, 231–273. Horne, R.N., Caltagirone, J.P., 1980. On the evaluation of thermal disturbances during natural convection in a porous medium. Journal of Fluid Mechanics 100, 385–395. Horton, C.W., Rogers, F.T., 1945. Convection currents in a porous medium. Journal of Applied Physics 16, 367–370. Kalia, N., Balakotaiahand, V., 2007. Modeling and analysis of wormhole formation in reactive dissolution of carbonate rocks. Chemical Engineering Science 62, 919–928. Kaviany, M., 1984. Thermal convective instabilities in a porous medium. ASME Journal of Heat Transfer 106, 137–142. Kuhn, M., Dobertb, F., Gessner, K., 2006. Numerical investigation of the effect of heterogeneous permeability distributions on free convection in the hydrothermal system at Mount Isa, Australia. Earth and Planetary Science Letters 244, 655–671. Lapwood, E.R., 1948. Convection of a fluid in a porous medium. Proceedings of the Cambridge Philosophical Society 44, 508–521. Leach, D.L., Rowan, E.L., 1986. Genetic link between Ouachita foldbelt tectonism and the Mississippi Valley-type lead-zinc deposits of the Ozarks. Geology 14, 931–935. Lebon, G., Cloot, A., 1986. A thermodynamical modeling of fluid flows through porous media: application to natural convection. International Journal of Heat and Mass Transfer 29, 381–390. Lencka, M.M., Anderko, A., Sanders, S.J., Young, R.D., 1997. Modeling viscosity of multicomponent electrolyte solutions. International Journal of Thermophysics 19, 367–378. Lichtner, P.C., 1996. Continuum formulation of multicomponent–multiphase reactive transport. In: Lichtner, P.C., Steefel, C.I., Oelkers, E.H. (Eds.), Reactive Transport in Porous Media, Mineralogical Society of America: Reviews in Mineralogy, 34, pp. 1–81. Lin, G., Zhao, C., Hobbs, B.E., Ord, A., Muhlhaus, H.B., 2003. Theoretical and numerical analyses of convective instability in porous media with temperature-dependent viscosity. Communications in Numerical Methods in Engineering 19, 787–799. Lister, J.R., 1990. Buoyancy-driven fluid fracture: the effects of material toughness and of low-viscosity precursors. Journal of Fluid Mechanics 210, 263–280. Lister, J.R., Kerr, R.C., 1991. Fluid-mechanical models of crack propagation and their application to magma transport in dykes. Journal of Geophysical Research 96, 10049–10077. Magri, F., Akar, T., Gemici, U., Pekdeger, A., 2010. Deep geothermal groundwater flow in the Seferihisar–Balçova area, Turkey: results from transient numerical simulations of coupled fluid flow and heat transport process. Geofluids 10, 388–405.

33

Marsh, B.D., 1982. On the mechanics of igneous diapirism, stoping, and zone melting. American Journal of Science 282, 808–855. Matthai, S.K., Henley, R.W., Heinrich, C.A., 1995. Gold precipitation by fluid mixing in bedding-parallel fractures near carbonaceous slates at the Cosmopolitan Howley gold deposit, northern Australia. Economic Geology and the Bulletin of the Society of Economic Geologists 90, 2123–2142. Matthai, S.K., 2003. Fluid flow and (reactive) transport in fractured and faulted rock. Journal of Geochemical Exploration 78–79, 179–182. Matthai, S.K., Geiger, S., Roberts, S.G., 2004. The Complex Systems Platform CSP5.0: User's Guide, ETH Research Reports, 5th edition. ETH Zurich, Zurich, Switzerland. Matthai, S.K., Nick, H.M., Pain, C., Neuweiler, I., 2010. Simulation of solute transport through fractured rock: a higher-order accurate finite-element finite-volume method permitting large time steps. Transport in Porous Media 83, 289–318. Miller, C.T., Christakos, G., Imhoff, P.T., McBride, J.F., Pedit, J.A., Trangenstein, J.A., 1998. Multiphase flow and transport modeling in heterogenous porous media: challenges and approaches. Advances in Water Resources 21 (2), 77–120. Mills, R.T., Lu, C., Lichtner, P.C., Hammond, G.E., 2007. Simulating subsurface flow and transport on ultrascale computers using PFLOTRAN. Journal of Physics Conference Series 78. doi:10.1088/1742-6596/78/1/012051. Mills, R.T., Lu, C., Hammond, G.E., Lichtner, P.C., Sripathi, V., Mahinthakumar, G., Smith, B.F., 2009. Modeling subsurface reactive flows using leadership-class computing. Journal of Physics Conference Series 180. doi:10.1088/1742-6596/180/1/012062. Molson, J.W., Frind, E.Q., 1993. Heatflow: Density-dependent Flow and Thermal Energy Transport Model in Three Dimensions, User Guide. Waterloo Center for Groundwater Research, University of Waterloo, Canada. Nield, D.A., Bejan, A., 1992. Convection in Porous Media. Springer-Verlag, New York. Norton, D., Knight, J., 1977. Transport phenomena in hydrothermal systems: cooling plutons. American Journal of Science 277, 937–981. Oliver, J., 1986. Fluids expelled tectonically from orogenic belts: their role in hydrocarbon migration and other geologic phenomena. Geology 14, 99–102. Oliver, N.H.S., 1995. Review and classification of structural controls on fluid flow during regional metamorphism. Journal of Metamorphic Geology 14, 477–492. Ord, A., Hobbs, B.E., Zhang, Y., Broadbent, G.C., Brown, M., Willetts, G., Sorjonen-Ward, P., Walshe, J., Zhao, C., 2002. Geodynamic modelling of the Century deposit, Mt Isa Province, Queensland. Australian Journal of Earth Sciences 49, 1011–1039. Owen, D.R.J., Feng, Y.T., de Souza Neto, E.A., Cottrell, M.G., Wang, F., Andrade Pires, F.M., Yu, J., 2004. The modelling of multi-fracturing solids and particular media. International Journal for Numerical Methods in Engineering 60, 317–339. Paluszny, A., Matthai, S.K., Hohmeyer, M., 2007. Hybrid finite element–finite volume discretization of complex geologic structures and a new simulation workflow demonstrated on fractured rocks. Geofluids 7, 186–208. Phillips, O.M., 1991. Flow and Reactions in Permeable Rocks. Cambridge University Press, Cambridge. Phillips, W.J., 1972. Hydraulic fracturing and mineralization. Journal of the Geological Society of London 128, 337–359. Pillatsis, G., Taslim, M.E., Narusawa, U., 1987. Thermal instability of a fluid-saturated porous medium bounded by thin fluid layers. ASME Journal of Heat Transfer 109, 677–682. Pokrovskii, V.A., 1999. Calculation of the standard partial molal thermodynamic properties and dissociation constants of aqueous HCl0 and HBr0 at temperatures to 1000 °C and pressures to 5 kbar. Geochimica et Cosmochimica Acta 63, 1107–1115. Pruess, K., 1992. TOUGH2 — a General-Purpose Numerical Simulator for Multiphase Fluid and Heat Flow, Report LBL-20700. Lawrence Berkeley Labortory, Berkeley, California. Raffensperger, J.P., Garven, G., 1995. The formation of unconformity-type uranium ore deposits: coupled hydrochemical modeling. American Journal of Science 295, 639–696. Rubin, A.M., 1995. Propagation of magma-filled cracks. Annual Review of Earth and Planet Science 23, 287–336. Saltzer, S.D., Pollard, D.D., 1992. Distinct element modeling of structures formed in sedimentary overburden by extensional reactivation of basement normal faults. Tectonics 11, 165–174. Schafer, D., Schafer, W., Kinzelbach, W., 1998. Simulation of reactive processes related to biodegradation in aquifers: 1, structure of the three-dimensional reactive transport model. Journal of Contaminant Hydrology 31, 167–186. Schaubs, P., Zhao, C., 2002. Numerical modelling of gold-deposit formation in the Bendigo– Ballarat zone, Victoria. Australian Journal of Earth Sciences 49, 1077–1096. Schmidt Mumm, A., Brugger, J., Zhao, C., Schacht, U., 2010. Fluids in geological processes: the present state and future outlook. Journal of Geochemical Exploration 106, 1–7. Schulenberg, T., Muller, U., 1987. An improved model for two-phase flow through beds of coarse particles. International Journal of Multiphase Flow 13 (1), 87–97. Scott, D.R., Stevenson, D.J., 1986. Magma ascent by porous flow. Journal of Geophysical Research 91, 9283–9296. Shock, E.L., Helgeson, H., 1988. Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: correlation algorithms for ionic species and equation of state predictions to 5 kb and 1000°C. Geochimica et Cosmochimica Acta 52, 2009–2036. Shock, E.L., Koretsky, C.M., 1995. Metal-organic complexes in geochemical processes: Estimation of standard partial molal thermodynamic properties of aqueous complexes between metal cations and monovalent organic acid ligands at high pressures and temperatures. Geochimica et Cosmochimica Acta 59, 1497–1532. Sibson, R.H., 1987. Earthquake rupturing as a hydrothermal mineralizing agent. Geology 15, 701–704. Sibson, R.H., Moore, J., Rankin, A.H., 1975. Seismic pumping: a hydrothermal fluid transport mechanism. Journal of the Geological Society 131, 653–659.

34

C. Zhao et al. / Journal of Geochemical Exploration 112 (2012) 21–34

Sibson, R.H., Robert, F., Poulsen, K.H., 1988. High angle reverse faults, fluid pressure cycling, and mesothermal gold-quartz deposits. Geology 16, 551–555. Sibson, R.H., 1994. Crustal stress, faulting and fluid flow. Geological Society Special Publication (No. 78), 69–84. Sorjonen-Ward, P., Zhang, Y., Zhao, C., 2002. Numerical modelling of orogenic processes and mineralization in the south eastern part of the Yilgarn Craton, Western Australia. Australian Journal of Earth Sciences 49, 935–964. Steefel, C.I., Lasaga, A.C., 1994. A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems. American Journal of Science 294, 529–592. Steefel, C.I., Yabusaki, S.B., 1996. OS3D/GIMRT: software for multicomponent– multidimensional reactive transport. User Manual and Programmer's Guide. PNL-11166. Pacific Northwest National Laboratory, Richland, Washington. Tournier, C., Genthon, P., Rabinowicz, M., 2000. The onset of natural convection in vertical fault planes: consequences for the thermal regime in crystalline basements and for heat recovery experiments. Geophysical Journal International 140, 500–508. Valenta, R.K., Cartwright, I., Oliver, N.H.S., 1994. Structurally controlled fluid flow associated with breccia vein formation. Journal of Metamorphic Geology 12, 197–206. Wang, C.Y., Beckermann, C., 1993. A two-phase mixture model of liquid–gas flow and heat transfer in capillary porous media: Formulation. International Journal of Heat and Mass Transfer 36 (11), 2747–2758. Wang, C.Y., Cheng, P., 1996. A multiphase mixture model for multicomponent transport in capillary porous media: model development. International Journal of Heat and Mass Transfer 39 (17), 3607–3618. Wang, P., Anderko, A., Young, R.D., 2002. A speciation-based model for mixed-solvent electrolyte systems. Fluid Phase Equilibria 203, 141–176. Wang, P., Anderko, A., 2008. Modeling thermal conductivity of concentrated and mixed-solvent electrolyte systems. Industrial and Engineering Chemistry Research 47, 5698–5709. Wang, P., Anderko, A., Springer, R.D., Kosinski, J.J., Lencka, M.M., 2010. Modeling chemical and phase equilibria in geochemical systems using a speciation-based model. Journal of Geochemical Exploration 106, 219–225. Weinberg, R.F., 1996. Ascent mechanism of felsic magmas: news and views. Transactions of the Royal Society of Edinburgh Earth Science 87, 95–103. Wilde, A.R., Wall, V.J., 1987. Geology of the Nabarlek Uranium deposit, Northernterritory, Australia. Economic Geology 82, 1152–1168. Wolery, T.J., 1979. Calculation of Chemical Equilibria between Aqueous Solution and Minerals: the EQ3/6 Software Package, UCR-52658. Lawrence Livermore Laboratory, Livermore, California, USA. Yang, J.W., 2006. Full 3D 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.W., Feng, Z., Luo, X., Chen, Y., 2010. 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, 236–243. Yang, J.W., 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.

Yeh, G.T., Tripathi, V.S., 1991. A model for simulating transport of reactive multispecies components: model development and demonstration. Water Resources Research 27, 3075–3094. Zhang, Y., Hobbs, B.E., Ord, A., Barnicoat, A., Zhao, C., Walsh, J.L., Lin, G., 2003. The influence of faulting on host-rock permeability, fluid flow and ore genesis of gold deposits: a theoretical 2D numerical model. Journal of Geochemical Exploration 78–79, 279–284. Zhang, Y., Roberts, P., Murphy, B., 2010. Understanding regional structural controls on mineralization at the Century deposit: a numerical modeling approach. Journal of Geochemical Exploration 106, 244–250. Zhao, C., Mühlhaus, 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., 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., Walshe, J.L., Mühlhaus, H.B., Ord, A., 2001. Finite element modelling of fluid–rock interaction problems in pore-fluid saturated hydrothermal/ sedimentary basins. Computer Methods in Applied Mechanics and Engineering 190, 2277–2293. Zhao, C., Hobbs, B.E., Ord, A., Peng, S., Mühlhaus, H.B., Liu, L., 2004. Theoretical investigation of convective instability in inclined and fluid-saturated three-dimensional fault zones. Tectonophysics 387, 47–64. Zhao, C., Hobbs, B.E., Hornby, P., Ord, A., Peng, S., 2006. Numerical modelling of fluids mixing, heat transfer and non-equilibrium redox chemical reactions in fluidsaturated porous rocks. International Journal for Numerical Methods in Engineering 66, 1061–1078. Zhao, C., Hobbs, B.E., Ord, A., Hornby, P., Peng, S., Liu, L., 2007a. Mineral precipitation associated with vertical fault zones: the interaction of solute advection, diffusion and chemical kinetics. Geofluids 7, 3–18. Zhao, C., Hobbs, B.E., Ord, A., Peng, S., Liu, L., 2007b. An upscale theory of particle simulation for two-dimensional quasi-static problems. International Journal for Numerical Methods in Engineering 72, 397–421. Zhao, C., Hobbs, B.E., Ord, A., 2008a. Convective and Advective Heat Transfer in Geological Systems. Springer, Berlin. Zhao, C., Hobbs, B.E., Ord, A., 2008b. Investigating dynamic mechanisms of geological phenomena using methodology of computational geosciences: an example of equal-distant mineralization in a fault. Science in China Series D: Earth Sciences 51, 947–954. Zhao, C., Hobbs, B.E., Ord, A., Peng, S., 2008c. Particle simulation of spontaneous crack generation associated with the laccolithic type of magma intrusion processes. International Journal for Numerical Methods in Engineering 75, 1172–1193. Zhao, C., Hobbs, B.E., Hornby, P., Ord, A., Peng, S., Liu, L., 2008d. Theoretical and numerical analyses of chemical-dissolution front instability in fluid-saturated porous rocks. International Journal for Numerical and Analytical Methods in Geomechanics 32, 1107–1130. Zhao, C., Hobbs, B.E., Ord, A., 2009. Fundamentals of Computational Geoscience: Numerical Methods and Algorithms. Springer, Berlin. Zhao, C., Hobbs, B.E., Ord, A., 2010. Theoretical and numerical investigation into roles of geofluid flow in ore forming systems: integrated mass conservation and generic model approach. Journal of Geochemical Exploration 106, 251–260.