Trans. Nonferrous Met. Soc. China 29(2019) 165-177

Numerical simulation of hydrothermal mineralization associated with simplified chemical reactions in Kaerqueka polymetallic deposit, Qinghai, China

Yan-hong ZOU^1,2, Yao LIU^1,3, Yong PAN⁴, Kuan-da YANG^1,2, Ta-gen DAI^1,2, Xian-cheng MAO^1,2, Jian-qing LAI^1,2, Hai-long TIAN⁵

1. Computational Geosciences Research Center, Central South University, Changsha 410083, China;

2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring, Ministry of Education, School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;

3. Embedded and Networking Computing Laboratory, Hunan University, Changsha 410082, China;

4. College of Resources Environment and Tourism, College of Hunan Arts and Science, Changde 415000, China;

5. Key Laboratory of Groundwater Resources and Environment, Ministry of Education, Jilin University, Changchun 130021, China

Received 19 November 2017; accepted 20 July 2018

Abstract: The Kaerqueka polymetallic deposit, Qinghai, China, is one of the typical skarn-type polymetallic ore deposits in the Qimantage metallogenic belt. The dynamic mechanism on the formation of the Kaerqueka polymetallic deposit is always an interesting topic of research. We used the finite difference method to model the mineralizing process of the chalcopyrite in this region with considering the field geological features, mineralogy and geochemistry. In particular, the modern mineralization theory was used to quantitatively estimate the related chemical reactions associated with the chalcopyrite formation in the Kaerqueka polymetallic deposit. The numerical results indicate that the hydrothermal fluid flow is a key controlling factor of mineralization in this area and the temperature gradient is the driving force of pore-fluid flow. The metallogenic temperature of chalcopyrite in the Kaerqueka polymetallic deposit is between 250 and 350 °C. The corresponding computational results have been verified by the field observations. It has been further demonstrated that the simulation results of coupled models in the field of emerging computational geosciences can enhance our understanding of the ore-forming processes in this area.

Key words: numerical simulation; heat transfer; chemical reaction; hydrothermal mineralization; Kaerqueka polymetallic deposit

1 Introduction

The Kaerqueka is located in Qimantage area of the western part of east Kunlun Mountain, Qinghai province, and the Kaerqueka polymetallic deposit is a large-scale characteristic hydrothermal vein-skarn ore deposit [1-3]. With the implementation of a new round of land and resources survey, the Kaerqueka polymetallic deposit has been studied extensively and intensively by many researchers [2,3]. In recent years, the studies [1-3] of the regional metallogenic regularity have made significant progress, but the related researches between the intermediate-acid intrusive magmatism and mineralization in this region are very limited, the mineralizing process of the chalcopyrite in Kaerqueka polymetallic deposit still needs further exploration. In this work, we used a new simulation method [4-6] to model the mineralizing process of the chalcopyrite and to get better understanding of ore-formation and mineralization in the Kaerqueka polymetallic deposit. The preliminary research results [1-3] provided some useful information about the Kaerqueka polymetallic deposit. We can briefly summarize the preliminary research results as follows: (1) the storage situation and distribution characters of the chalcopyrite resources in the Kaerqueka polymetallic deposit have been investigated; (2) the strata and the structure in this district have been determined; (3) the metallogenetic resource comes from the hydrothermal fluids; (4) the common way of the ore forming is through the hydrothermal fluids to penetrate a carbonate layer, where a chemical reaction should take place.

Hydrothermal mineralization is an extremely complicated process, and the traditional metallogenic theories for revealing the ore-forming process have some limitations [7-9]. Because there are many factors that can affect the formation of mineral deposit, such as the pore-fluid flow channel, temperature distribution, pore pressure variation and physical chemical properties of surrounding rocks [10-12], it is hard for the traditional metallogenic theories to obtain sufficient details associated with an ore-forming system as well as the critical geodynamic processes responsible for formation and localization of orebodies [13-15]. Therefore, the numerical simulation method becomes a powerful tool for overcoming their drawbacks [16-18]. With the significant advancements of computational science and computer technology, the computational simulation method, as a new method for dealing with a lot of complex problems in the field of mineral exploration, has attracted increasing interest of research. Extensive studies of the ore-forming processes in the upper crust of the earth have been conducted by ZHAO et al [19-21] and others [22-24]. The previous studies can be briefly classified into the following categories: (1) numerical simulation and theoretical analyses of coupled heat transfer and pore-fluid flow processes without considering the rock deformation [25-27]; (2) numerical simulation and theoretical analyses of convective and advective heat transfer in geological systems [28-30];  (3) theoretical study of chemical dissolution instability [31-33] and mineral precipitation problems [34-36]; (4) numerical simulation of feedback effects of chemical reactions on porosity and fluid flow channels [37-39]; and (5) application of numerical simulation to different ore-forming systems in different geological regions [40-42] and also to various types of geoscience problems [43-45].

However, the ore-forming process in the Kaerqueka polymetallic deposit is very complex. There are still several key factors that have been ignored. For example: (1) the effects of temperature and pressure gradients on the chalcopyrite ore deposits have not been analyzed in this district; (2) the processes of ore migration are seldom considered; (3) the researches of related chemical reactions for chalcopyrite ore precipitation are still not deep enough; (4) many studies are restricted to the traditional methods, and the new approach of computational simulation [4,17] has not been tried in the Kaerqueka polymetallic deposit; and (5) the chemical factor of the chalcopyrite ore has seldom been analyzed. Therefore, it is necessary to develop innovative knowledge models for understanding ore-forming processes in detail, so as to facilitate the under-standing of ore-formation processes in the Kaerqueka polymetallic deposit.

Among above-mentioned research  achievements, ZHAO et al [46-48] have conducted primitive and pioneering work on reactive fluid-mixing and mineralization in pore-fluid hydrothermal systems. In particular, ZHAO et al [12,17,49] proposed the mineralization rate concept and modern mineralization theory, which were successfully applied to the finite element simulation of mineralization patterns in hydrothermal systems, and excellent results have been obtained [50-52]. In this study, we used the numerical simulation method [4,53] and the modern mineralization theory [49] to investigate the chalcopyrite ore-forming mechanism of the Kaerqueka polymetallic deposit. The outcome of this study will not only enrich the contents of the emerging computational geoscience discipline [54-56], the methodology of which has been used to solve many geoscience problems [57-60], but also give us a better understanding of the ore-formation processes in this particular region.

2 Geological setting and metallogetic model of ore deposit

2.1 Geological setting

The Qimantage region in Qinghai Province is located in the southwest margin of the Qaidam Basin within the northeast of Qinghai-Tibet Plateau and lies in the Qaidam block along the east Kunlun orogenic belt, which is known as an important polymetallic metallogenic belt [1]. The Kaerqueka polymetallic deposit, located in the Qimantage area of the western part of the east Kunlun orogenic belt, is an important polymetallic ore deposit area in Qinghai Province. It is also the main mining area for the project of the integrated exploration in Kaerqueka [61]. The deposit was discovered by Qinghai Geological Survey Institute, China in 2003. The ore-forming elements are complicated, which are dominated by Cu, Pb and Zn and accompanied with Fe, Au, Ag and so on [61]. The exposed rocks were mainly granitic intrusive rocks, followed by a small amount of Ordovician-Silurian Tanjianshan group of marble and basic volcano rocks [2]. In recent years, the amount of copper polymetallic resources in mining areas has reached 5×10⁴ t, and has become the large-scale deposit. The mineralization was dominated by skarn-type polymetallic ore and had a feature of porphyry-skarn-hydrothermal series [3].

2.1.1 Strata

The strata from early to late time sequences in this area can be found in the previous reference [62]. The carbonate rocks were present in all strata, especially developed in the Dagangou and Di’aoshu Formation in carboniferous. The tectonic activities in the region were developed, mainly dominated by the NW fault structure and NWW fault structure, with the NE and nearly NS fault structure formed the main body of the tectonic framework, which controlled the distribution of strata, magmatic rocks and minerals [62]. Magmatic activities were frequent and dominated by Hercynian-Indosinian magmatism, especially strong in the intermediate acid magma intrusion. The magmatic activities were closely related to the mineralization of the porphyry and skarn deposit in this area and provided material basis and energy for mineralization [62].

The exposed strata in this mining area is relatively simple, including the Palaeoproterozoic Jinshuikou group Baishahe rock group, Ordovician-Silurian Tanjianshan group, Elashan Formation of upper Triassic series and quaternary [3]. Palaeoproterozoic Baishahe Formation of the Jinshuikou group was mainly distributed in the northern part of the mining area, and its lithology was mainly composed of biotite monzonite gneiss. The Ordovician-Silurian Tanjianshan group was surrounded by rock mass, which formed an island-like remnant [63]. Its lower part was the basic volcanic rock primarily composed of andesite-basalt and amphibolite, and the upper part is the combination of dolomite, marble, slate and phyllite. Elashan Formation of upper Triassic series was sporadically distributed in the mouth of the valley and the front of the hillside of Kaerqueka, where the angle was not integrated with the stratum of Palaeoproterozoic Jinshuikou group [3], and its lithology was made up of rhyolite-tuff, rhyolitic breccia ignimbrite, andesite, and volcanic breccia [63]. Quaternary was chiefly combined by moraine, alluvial gravel, alluvial gravel, silty sand and sub clay [1].

2.1.2 Structure

The structures of the mining area were mainly NWW and NW faults. The former was the main ore-controlling structure of the mining area [63], and developed in the Late Permian epoch’s biotite monzogranite, which was a group of near-parallel faults with compression, torsion, multi-stage activities and so on. The NE faults were small in size and developed in the late Triassic granodiorite, which cut the NW faults and the Late Permian granitic batholith.

2.1.3 Magmatic rock

The magmatic rocks in the mining area were complex, intermediate-acidic and basic-ultrabasic rocks were exposed [2], the intrusive rocks were dominated  by porphyroid biotite adamellite and granodiorite. Porphyroid biotite adamellite occurred as the batholith and strike along the NNW direction in the whole, which was basically consistent with regional tectonic lines. Granodiorite intruded into the former in the form of rock plant, and the shape was regular long strip and was also controlled by NWW with obvious structural control. In addition, some diorite, diorite-porphyrite and granite were also distributed as small rocky branch or vein [63].

2.1.4 Mineralization

The mineralization was complex and diverse. There were 64 orebodies about copper, zinc and silver. According to the mineralization type, the ore can be divided into A, B, C and D zones (Fig. 1). Among them, zones A, B and C are the main mineralized areas [64].

Fig. 1 Geological sketch map of Kaerqueka polymetallic polymetallic ore district [63]

The porphyry mineralization mainly occurred in three nearly parallel cataclastic alteration zones in zone A, which extended in the NNW direction. The main lithology in the cataclastic alteration zone was porphyroid biotite adamellite, which was filled by later fine grained granite and quartz vein. The skarn mineralization belt was mainly located in the contact area between granodiorite and Tanjianshan group. Part of the porphyritic granite contact zone in the Late Permian was also mineralized and with zonal distribution. The interior of the zonal distribution was composed of diopside skarn, tremolite skarn, diopside garnet, skarn, marble, and andesite. The alteration mainly consists of skarnization, chloritization, epidotization, silicification and so on and the main ore minerals were chalcopyrite, chalcocite, bornite, molybdenite, magnetite, hematite and so on [64].

2.2 Metallogetic model for Kaerqueka

The metallogenic belt was formed in the eastern Kunlun collisional orogenic belt during the magmatic collision and underplating in Late Indosinian epoch and produced in the intersecting parts between the NWW trending faults associated with the middle and deep faults of the Kunlun Mountains and NE trending fracture in Altyn Tagh [63].

In the late Indosinian epoch, due to the subduction of the south Kunlun landmass to middle Kunlun landmass, the large-scale inland orogeny was formed, the crust was thickened and gradually changed from the main collision to the backward collision, and the stress mechanism was changed from extrusion to extension. Subduction caused partial melting of the lower crust to form magma. The stress background of the crustal thickening and relative relaxation caused the delamination of lower crust, upwelling of the mantle material, and contamination of magma fused with the crust. The blazing magma emplaced upwards along the favorable transport path and carried a large number of mineral elements. At the same time, the metal elements containing in different lithologic sections of the strata of the Tanjianshan group are different. During the process of intrusion, the magma underwent material metasomatism with the crust of Tanjianshan group and formed skarn deposits near or between the contact zones. Due to different lithologies of metasomatic surrounding rock, the lithology of zone B was volcanic and marble, and the copper-iron-zinc mineralization was dominated, while the surrounding rock of zone C was marble and the lead-zinc-silver mineralization was dominated. During the ascent process of magma, owing to the change of physical and chemical conditions, the ore-forming elements were enriched in the rock mass to form a porphyry deposit. The ore-forming fluids moved along the crushing zone, enriched and precipitated in the favorable parts of the crushing zone, that is the reason for the formation of a shallow-type low-temperature hydrothermal deposit [65].

The metallogenic model of Kaerqueka polymetallic deposit can be expressed by the following diagram (Fig. 2 [66]).

Fig. 2 Metallogenic model of Kaerqueka polymetallic deposit [66]

3 Geological/conceptual models for Kaer- queka polymetallic deposit

Fig. 3 Generic model of Kaerqueka copper polymetallic deposit

We constructed two computational models, a generic model and a realistic model, for simulating the thermal flow driven mineralization processes associated with the Kaerqueka polymetallic deposit. A conceptual 2D model (Fig. 3) is established for the Kaerqueka polymetallic deposit based on the cross-section along line EF in Fig. 1. A set of assumptions are made to reflect the subsurface geological conditions of the Kaerqueka polymetallic deposit and other features that are relevant to pore-fluid flow in this district. In this two models, the depth of the model is 2 km; the horizontal dimension of the model is 1 km; the driving force of pore fluid flow is the heat anomaly and the time of simulation is set to be 10000 years; the top surface is set to be 25 °C and kept to be constant. The initial hot pore-fluid in the fault zone is 650 °C, because the metallogenic temperature is within 500 °C, and there is basically no ore-forming chemical reaction over 500 °C in this region. In addition, before the ore-forming chemical reaction takes place, the ore-bearing hydrothermal fluid needs to flow a certain distance. For this reason, we deliberately raised the initial temperature to 650 °C in the models. The pore fluid pressure of both the sedimentary rocks and skarn is set to be hydrostatic pressure, while the initial pore fluid pressure of the intrusion is set to be lithostatic pressure. These chosen conditions are based on the following considerations. The existing study [62] showed that the deepest depth of the exploration line is almost 1 km. The particular distance can be calculated by the depth of the exploration line with considering the distance of fluid flow. We can assume that the depth of the model is 2 km. The horizontal dimension of the model is 1 km, which is similar to the length of line CD in Fig. 1. As mentioned previously, the mineralization temperature of the chalcopyrite is 225-375 °C. Since the highest temperature used for the TOUGHREACT calculation in the experimental data of the equilibrium concentration of copper ion (Cu²⁺) is about 400 °C, the equilibrium concentration of copper ion at temperatures between 400 to 500 °C is the extrapolated data. The initial temperature of the intrusion is set to be 650 °C because the initial temperature of the intrusion in the experiments must be a little higher than the highest temperature of ore-forming fluids [7]. Thus, 650 °C is an acceptable initial temperature for the intrusion in the following numerical simulations. As above mentioned, the structure was formed before the mineralization, and the deformation did not participate in the ore-forming process [62]. When the intrusion arrived at the bottom of the fault, the temperature difference provided the driving force of the pore-fluid flow, so we believed that the driving force of the pore-fluid flow is the heat anomaly.

The 2D model consists of 4800 four-node grid cells. According to the available geological data [62], many orebodies are located in the contact zone between the magmatic intrusion and the carbonate layer. As shown in Fig. 3, the main wall-rock types are made up of marble, skarn, quaternary, granodiorite and fault. The depth of the fault zone is assumed to be almost the perpendicular distance. Because many orebodies are exposed on the surface [3], which implies that the fault may extend to the surface, the depth of fault zone was determined based on this recognition. Since this paper was only interested in the temperature dependence of mineralization, the rock deformation of the intrusion was neglected. The relevant rock material parameters in Table 1 are used in the numerical simulation.

The fault distributions in the model are based on the geological cross-section map in Fig. 4. To consider the rock erosion situation during the mineralization and the fluid mobility, the extended model is established and filled with the corresponding rock. The numerical model of the second example is constructed by considering one of the geological section maps of the Kaerqueka polymetallic deposit, as shown in Fig. 5. This model, which is called the “realistic” geological conceptual model, is filled with more complex rocks and composed of 9600 four-node grid cells. In Fig. 5, there are two parallel fault zones and the dip angles both being 65° to the east.

Table 1 Rock material parameters assigned for generic model (Data are obtained from previous studies [4,43,53])

Fig. 4 Geological section along zero exploration line in mine district B of Kaerqueka deposit [2]

Fig. 5 “Realistic” geological conceptual model of Kaerqueka polymetallic deposit [2]

4 Establishment of mathematical model for Kaerqueka polymetallic deposit

4.1 Brief introduction of FLAC software

The Fast Lagrangian Analysis of Continua (FLAC) code [4,67] is used to simulate the mineralization process of the Kaerqueka polymetallic deposit. FLAC2D is a two-dimensional explicit finite difference code, which has been widely used for solving complex geoscience problems [68-70].

Compared with the finite element models, the explicit finite difference models have the following advantages: (1) it can be used to calculate nonlinear constitutive relationships directly, without excessive memory requirements; (2) a large number of grid cells may be modeled with a modest memory requirement, because there is no global stiffness matrix that needs to be inverted during each time-step; (3) FLAC2D is robust in the sense and can handle any constitutive model with no adjustment to the solution algorithms and techniques for different constitutive models. The open programming design option makes it easy for the users to use their own programming techniques within the framework of the FLAC2D code through the FISH language. Large-scale complex problems can be simulated using the FLAC2D code, and computational time would not increase significantly for highly nonlinear problems [67,68].

4.2 Governing mathematical equations associated with mathematical model for Kaerqueka copper polymetallic deposit

As mentioned previously, the hydrothermal trans- formation of the chalcopyrite should be caused by a thermal process, which is controlled by the temperature gradient in porous rocks [12,17]. As the hydrothermal fluid passed through the carbonate layer, a strongly chemical reaction occurred, so that the chalcopyrite was formed. Since the mineralization rate, which is presented by ZHAO et al [12,49], has been successfully used to simulate the precipitation of the magnetite in the Fushan ore deposit [4], it can be further extended to explain the other metals ore-forming process, in which the chalcopyrite was formed by hydrothermal transformation in the Kaerqueka polymetallic deposit. The governing equations used for simulating the chalcopyrite ore-forming process of the Kaerqueka polymetallic deposit in the FLAC models can be briefly described as follows:

                            (1)

                     (2)

                    (3)

                    (4)

                        (5)

where  and  are the Darcy velocities of fluid flow in x and y directions, respectively; P is the pressure of the pore-fluid; k_x and k_y are the mobility coefficients (permeability of FLAC) in x and y directions, respectively; λ(s) is a function of saturation s; μ is the dynamic viscosity of the pore-fluid; g is the component of gravitational acceleration in y direction; ρ_f is the density of the pore-fluid;  and  are the thermal flux values in x and y directions, respectively;  and  are the thermal conductivity coefficients of the pore-fluid and solid matrix, respectively;  is the porosity; ρ₀ is the reference density of the pore-fluid; β is the volumetric thermal expansion coefficient; T and T₀ are the temperature and reference temperature of the pore-fluid and surrounding medium, respectively.

Equations (1) and (2) indicate the Darcy law describing pore-fluid flow in x and y directions, respectively; Eqs. (3) and (4) indicate the Fourier law describing heat transfer in x and y directions, respectively; while Eq. (5) shows the relationship between temperature and the pore-fluid density [71-73]. The governing equations describe the conservation of mass, energy and momentum. The coupled thermal-hydrological constitutive relations can be found in Refs. [74-76].

4.3 Concept of mineralization rate

The mineralization rate of a mineral, which was proposed by ZHAO et al [12,49], is defined as the variation in the mineral mass per unit volume rock at per unit time during mineralization. Using this definition, the positive value of the mineralization rate of copper ore means the dissolution of the chalcopyrite in the hydrothermal system, while the negative value of the mineralization rate means the precipitation of the copper and iron ores in the hydrothermal system. The mathematical formula of the mineralization rate can be deduced from the mass conservation law as   follows [12,49]:

(k=1, 2, …, M)          (6)

where u_j is the velocity component in the x_j direction;  is the equilibrium concentration of mineral k; R_k is the source/sink term; M is the total number of the minerals to be considered in the ore-forming system; andis the second-order diffusivity of the porous rock. For mineralization problems, the diffusion term on the right-hand side of Eq. (6) is usually much smaller than the advection term. Thus, Eq. (6) can be approximately expressed as

                          (7)

where is the mineralization rate of mineral k.

In general, the equilibrium concentration of a mineral is a function of temperature, pressure and other relevant minerals as follows:

=f(T, P, C₁, C₂, …, C_N)                     (8)

where C_r (r=1, 2, …, N) is the concentration of mineral; N is the number of the relevant minerals/species to determine the equilibrium concentration of mineral/ species k in the chemical reaction. According to the chain rule, substituting Eq. (8) into Eq. (7) yields the following equation [12,49]:

  (9)

4.4 Equilibrium concentration of copper ion

According to the previous studies [2,3], the chalcopyrite in the Kaerqueka ore district is mainly located in the contact belt between intermediate-acidic intrusive rocks and dolomite-marble. The main orebody is composed of chalcopyrite in Kaerqueka ore district. Since the equilibrium concentration of a mineral is strongly dependent on the local temperature, the temperature variation in a hydrothermal system may provide a favorable environment for the mineralization to occur.

For the mixing of reactive sulfide fluids, the corresponding chemical formula can be expressed as

                          (10)

          (11)

                  (12)

where w is the overall reaction rate constant, which    is strongly dependent on the temperature, from the thermodynamic point of view. The chemical formulae [Cu(NH₃)₄](OH)₂ and Na[Fe³⁺Cl₄] are the common complexes of copper ion in the hydrothermal fluids.

The previous studies [2,77,78] showed that the chalcopyrite could be easily precipitated when the fluid is neutral solution. Because the pH value of NaCl aqueous solution almost equals 7, it can be used as a buffer for the chalcopyrite mineralization. Since the equilibrium concentration of copper ion in chalcopyrite has been directly calculated by the TOUGHREACT  code [79] in NaCl solution that comprises of calcium carbonate, sodium chloride, the complex of copper and iron, hydrogen sulfide, and ferrous chloride, it can be used in this study for analyzing the precipitation pattern of the chalcopyrite. For a NaCl buffer, the equilibrium concentration of copper ion in chalcopyrite between 25 and 500 °C can be expressed in a piecewise manner as

  (13)

where  is the equilibrium concentration of copper ion. Note that the equilibrium concentrations of copper ion at temperatures between 400 and 500 °C are the extrapolated data.

Figure 6 shows the variation of copper ion equilibrium concentration with temperature under the NaCl buffer solution condition, while Fig. 7 shows the variation of copper ion equilibrium concentration with pressure under the NaCl buffer solution condition. Similarly, Fig. 8 shows the variation of the differentiation of copper ion equilibrium concentration with temperature under the NaCl buffer solution condition. The solid line is the experimental data, while the dotted line is the extrapolated data.

Fig. 6 Variation of copper ion equilibrium concentration with temperature in NaCl buffer (Experimental data are provided by TOUGHREACT code [79])

Fig. 7 Variation of copper ion equilibrium concentration with pressure in NaCl buffer solution (Experimental data are provided by TOUGHREACT code [79])

Fig. 8 Variation of differentiation of copper ion equilibrium concentration with temperature in NaCl buffer solution (Experimental data are provided by TOUGHREACT code [79])

Based on Figs. 6-8, the experimental data of copper ion concentration were obtained under the conditions that the concentrations of sulfide ion, sulfate ion and carbonate ion are not constant. By using the TOUGHREACT-calculated equilibrium concentrations of copper ion, three-valent iron ion, sulfide, sulfate ion, and carbonate ion at different temperatures and  pressures, it is obvious that the equilibrium concentration of copper ion, three-valence iron ion, sulfide, sulfate ion, and carbonate ion are controlled by temperature rather than by pressure, because pressure has little effect on the equilibrium concentration of copper ion.

5 Related simulation results for Kaerqueka copper polymetallic deposit

Since this work is only interested in the temperature dependence of mineralization, and the structures have been formed before the mineralization, the rock deformation of the intrusion is neglected. The relevant rock material parameters in Table 1 are used in the numerical simulation, and the simulation results of generic model are shown below.

Figure 9 shows the distributions of the temperature and the first-order partial derivative of the temperature with respect to the x and y directions in the generic model. As mentioned by ZHAO et al [12,17], the localization of temperature can provide appropriate conditions for the orebody formation in some regions of the hydrothermal system. From the mineralization rate theory, the place where more changes of temperature and the faster velocity of pore-fluid flow are more favorable to form orebodies [12,49]. Since the first-order partial derivative of temperature with respect to x and y directions represents the temperature gradients in x and y directions, we can use them to roughly estimate the favorable locations to form the orebodies in the hydrothermal system.

Fig. 9 Distributions of temperature (a), temperature gradients in x (b) and y (c) directions of generic model

Figure 10(a) shows the approximate mineralization rate, while Fig. 10(b) shows the modified mineralization rate. The distributions of the approximate mineralization rate (which is equal to the scalar product of the pore-fluid velocity and temperature gradient [12,49] and the modified mineralization rate (which is equal to the product of the approximate mineralization rate and the equilibrium concentration differentiation of copper ion with respect to temperature [12] in the generic model. This clearly shows that the approximate mineralization rate has not considered the equilibrium concentration differentiation of copper ion, but the modified mineralization rate has taken into account the equilibrium concentration differentiation of copper ion. This is the main difference between them. Since the negative value of the mineralization rate of a mineral reflects the precipitation of the mineral and the positive value of the mineralization rate reflects the dissolution of the mineral, the precipitation regions of the chalcopyrite are clearly shown in the generic model. As can be seen from these modeling results, most precipitates of copper ore are distributed in the region between 250 and 350 °C. As observed in the Kaerqueka polymetallic deposit. Since the pressure dependence of solubility is typically weak and usually neglected in the chemical reaction of chalcopyrite, the change of pressure has not been coded in the program [4]. Because the experimental data of the equilibrium concentration of copper ion, which is provided by the TOUGHREACT code, shows that the pressure has also little impact to the equilibrium of the chemical reaction of the chalcopyrite, the influence of pressure on the chalcopyrite precipitation has not been considered in the related numerical simulations.

Fig. 10 Distributions of approximate (a) and modified (b) mineralization rates with only considering effect of temperature

The simulation results of geological conceptual model are shown below, and the relevant rock material parameters listed in Table 1 are used in the numerical modeling, also without considering the rock deformation.

Fig. 11 Distributions of temperature (a), temperature gradients in x (b) and y (c) directions of “realistic” geological model

Figure 11 shows the distributions of temperature and temperature gradient in x and y directions of the “realistic” geological model, respectively. The results in Fig. 9 show that the distribution of temperature is almost symmetrical in the fault trend direction completely, while Fig. 11 shows that because the location of the fault is near the west of the bottom, and the attitudes of the two faults have a dip angle of 65° to the east, the distribution of temperature is diffused along the fault trend. The temperature gradient in x and y directions shows that the temperature has the greatest variation in the regions which are favorable for the chalcopyrite precipitation. This indicates that both the position and intrusion directions of the fault zone have a significant influence on the temperature localization in this hydrothermal system.

Figure 12 shows the distributions of the approximate and modified mineralization rates of the chalcopyrite in the “realistic” geological model of the Kaerqueka copper polymetallic deposit. Figure 12(a) shows the approximate mineralization rate, while Fig. 12(b) shows the modified mineralization rate. Generally, the approximate mineralization rate can be only used to qualitatively predict the minaralization patern in the Kaerqueka polymetallic ore district, but cannot be used to quantitatively predict the specific amount of precipitated ore in the Kaerqueka copper polymetallic ore district. This is because the approximate mineralization rate is independent of a specific mineral [12,49]. On the contrary, since the modified mineralization rate is dependent on a specific mineral, it can be used to predict the specific amount of precipitated ore, if the mineralization time in the Kaerqueka polymetallic ore district is known.

Fig. 12 Distributions of approximate (a) and modified (b) mineralization rates of “realistic” geological model

6 Conclusions

(1) Hydrothermal fluid flow is one of the key controlling factor of mineralization in the Kaerqueka polymetallic deposit, and the driving force of pore-fluid flow is the temperature gradient. The metallogenic of the chalcopyrite mainly focuses on the big variation area of temperature and the region of pore-fluid flow with the faster velocity.

(2) According to the experimental data of the chemical reaction of the chalcopyrite, the pore-fluid pressure has little effect on the precipitation of chalcopyrite. The equilibrium concentration of copper ion has changed a lot between 250 and 350 °C. The ranges of metallogenic temperature of chalcopyrite, which are provided by the TOUGHREACT code, basically meet the field observations.

(3) Since the numerical results can reflect the field observations in this district basically, it has been further demonstrated that both the coupled thermao-hydrological model of the Kaerqueka polymetallic deposit and the numerical method for simulating pore-fluid flow and heat transfer associated with the chalcopyrite in this district are reasonable and feasible, although the simplified chemical reactions are considered in this research.

References

[1] LI Dong-sheng, GU Feng-bao, ZHANG Hai-lan, SU Sheng-shun, ZHONG Liang-yan, LIU Guang-liang. Geologic characteristics of the Kaerqueka porphyry copper deposit in Qinghai Province and its prospecting significance [J]. Northwestern Geology, 2012, 45: 174-183. (in Chinese)

[2] YU M, FENG C Y, ZHAO Y M, LI D X, XIAO Y, LIU J N, LI Z F. Fluid inclusion geochemistry in the Kaerqueka copper polymetallic deposit, Qinghai province and its genetic significances [J]. Acta Geologica Sinica, 2014, 88: 903-917. (in Chinese)

[3] GAO Y B, LI K, QIAN B, LI W Y, LI D S, SU S, ZHANG C G. The genesis of granodiorites and dark enclaves from the Kaerqueka deposit in east Kunlun belt: Evidence from zircon U-Pb dating, geochemistry and Sr-Nd-Hf isotopic compositions [J]. Geology in China, 2015, 42: 646-662. (in Chinese)

[4] ZOU Y H, LIU Y, DAI T G, MAO X C, LEI Y B, LAI J Q, TIAN H L. Finite difference modeling of metallogenic processes in the Hutouya Pb-Zn deposit, Qinghai, China: Implications for hydrothermal mineralization [J]. Ore Geology Reviews, 2017, 91: 463-476.

[5] SCHAUBS P, HOBBS B, ZHAO C. Computational simulation of seepage instability problems in fluid-saturated porous rocks: Potential dynamic mechanisms for controlling mineralization patterns [J]. Ore Geology Reviews, 2016, 79: 180-188.

[6] SCHAUBS P, HOBBS B, ZHAO C. Acquisition of spatially-distributed geochemical data in geoinformatics: Computational simulation approach [J]. Journal of Geochemical Exploration, 2016, 164: 18-27.

[7] LIU L, WAN C, ZHAO Y. Geodynamic constraints on orebody localization in the Anqing orefield, China: Computational modeling and facilitating predictive exploration of deep deposits [J]. Ore Geology Reviews, 2011, 43: 249-263.

[8] LUO Hai-jun, JIE Wan-qi, GAO Zhi-ming, ZHENG Yong-jian. Numerical simulation for macrosegregation in direct-chill casting of 2024 aluminum alloy with an extended continuum mixture model [J]. Transactions of Nonferrous Metals Society of China, 2018, 28: 1007-1015.

[9] SCHAUBS P, HOBBS B, ZHAO C. Effects of porosity heterogeneity on chemical dissolution-front instability in fluid-saturated rocks [J]. Journal of Central South University, 2017, 24: 720-725.

[10] ORD A, HOBBS E B, ZHAO C. Chemical dissolution-front instability associated with water-rock reactions in groundwater hydrology: Analyses of porosity-permeability relationship effects [J]. Journal of Hydrology, 2016, 540: 1078-1087.

[11] ZHAO C, HOBBS E B, ORD A. Modeling of mountain topography effects on hydrothermal Pb-Zn mineralization patterns: Generic model approach [J]. Journal of Geochemical Exploration, 2018, 190: 400-410.

[12] ZHAO C, HOBBS B E, ORD A. Convective and advective heat transfer in geological systems [M]. Berlin: Springer, 2008.

[13] ORD A, HOBBS E B, ZHAO C. Effects of acid dissolution capacity on the propagation of an acid-dissolution front in carbonate rocks [J]. Computers and Geosciences, 2017, 102: 109-115.

[14] HOBBS E B, ORD A, ZHAO C. Why asymptotic limit of the acid dissolution capacity can lead to a sharp dissolution front in chemical dissolution of porous rocks [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2017, 41: 1590-1602.

[15] ZHAO C, HOBBS B E, MUHLHAUS H B. Finite element modelling of temperature gradient driven rock alteration and mineralization in porous rock masses [J]. Computer Methods in Applied Mechanics and Engineering, 1998, 165: 175-187.

[16] HOBBS B E, ORD A. Investigating dynamic mechanisms of geological phenomena using methodology of computational geosciences: An example of equal-distant mineralization in a fault [J]. Science in China Series D: Earth Sciences, 2008, 51: 947-954.

[17] ZHAO C, HOBBS B E, ORD A. Fundamentals of computational geoscience: Numerical methods and algorithms [M]. Berlin: Springer, 2009.

[18] LIN G, HOBBS B E, ORD A, MUHLHAUS H B. Theoretical and numerical analyses of convective instability in porous media with temperature-dependent viscosity [J]. Communications in Numerical Methods in Engineering, 2003, 19: 787-799.

[19] ZHAO C, HOBBS E B, ORD A. Analytical solution for dissolution-time-scale reactive transport in fluid-saturated porous rocks [J]. International Journal of Geomechanics, 2018, 18: 1-12.

[20] ZHAO C, HOBBS E B, ORD A. Validity of using large-density asymptotics for studying reaction-infiltration instability in fluid-saturated rocks [J]. Journal of Hydrology, 2018, 559: 454-460.

[21] ZHAO C, HOBBS B, ORD A. A unified theory for sharp dissolution front propagation in chemical dissolution of fluid-saturated porous rocks [J]. Science China: Technological Sciences, 2018, 61: 1-12.

[22] LIN G, PENG M, ZHANG L, ZHANG D. Numerical analysis and simulation experiment of lithospheric thermal structures in the South China Sea and the Western Pacific [J]. Journal of Earth Science, 2009, 20: 85-94.

[23] ORD A, HOBBS B E, ZHANG Y, BROADBENT G C, BROWN M, WILLETTS G, SORJONEN-WARD P. Geodynamics modeling of the Century deposit, Mt Isa Province, Queensland, Australia [J]. Australian Journal of Earth Sciences, 2002, 49: 1011-1039.

[24] LIU L M, YANG G Y, PENG S L. Numerical modeling of coupled geodynamical processes and its role in facilitating predictive ore discovery: An example from Tongling, China [J]. Resource Geology, 2005, 55: 21-31.

[25] KUHN M, MUHLHAUS H B, PENG S. Numerical simulation of double-diffusion driven convective flow and rock alteration in three-dimensional fluid-saturated geological fault zones [J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195: 2816-2840.

[26] HORNBY P, PENG S, LIU L M. Mineral precipitation associated with vertical fault zones: The interaction of solute advection, diffusion and chemical kinetics [J]. Geofluids, 2007, 7: 3-18.

[27] HOBBS B E, ZHANG Y, ORD A. Application of coupled deformation, fluid flow, thermal and chemical modelling to predictive mineral exploration [J]. Journal of Geochemical Exploration, 2000, 70: 505-509.

[28] HOBBS B E, ORD A, PENG S, MUHLHAUS H B, LIU L M. Theoretical investigation of convective instability in inclined and fluid-saturated three-dimensional fault zones [J]. Tectonophysics, 2004, 387: 47-64.

[29] GOW P A, UPTON P, HILL K C. Copper-gold mineralisation in New Guinea: Numerical modelling of collision, fluid flow and intrusion-related hydrothermal systems [J]. Australian Journal of Earth Sciences, 2002, 49: 753-771.

[30] ORD A, PENG S, LIU L. Inversely-mapped analytical solutions for flow patterns around and within inclined elliptic inclusions in fluid-saturated rocks [J]. Mathematical Geosciences, 2008, 40: 179-197.

[31] ORD A, HORNBY P, PENG S. Morphological evolution of three-dimensional chemical dissolution front in fluid-saturated porous media: A numerical simulation approach [J]. Geofluids, 2008, 8: 113-127.

[32] ZHAO C B. Advances in numerical algorithms and methods in computational geosciences with modeling characteristics of multiple physical and chemical processes [J]. Science China: Technological Sciences, 2015, 58: 783-795.

[33] ORD A, HOBBS B E, ZHAO B C. Theoretical and numerical investigation into roles of geofluid flow in ore forming systems: Integrated mass conservation and generic model approach [J]. Journal of Geochemical Exploration, 2010, 106: 251-260.

[34] ORD A, HOBBS B E, ZHAO B C. Effects of medium and pore-fluid compressibility on chemical-dissolution front instability in fluid-saturated porous media [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2012, 36: 1077-1100.

[35] ZHAO C, HOBBS B E, ORD A. Theoretical analyses of acidization dissolution front instability in fluid-saturated carbonate rocks [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2013, 37: 2084-2105.

[36] HORNBY P, ORD A, PENG S, LIU L. Theoretical and numerical analyses of chemical-dissolution front instability in fluid-saturated porous rocks [J]. International Journal for Numerical and Analytical Methods in Geomechanics 2008, 32: 1107-1130.

[37] HOBBS B E, ORD A, HORNBY P, PENG S. Effect of reactive surface areas associated with different particle shapes on chemical-dissolution front instability in fluid-saturated porous rocks [J]. Transport in Porous Media, 2008, 73: 75-94.

[38] HOBBS B E, ORD A, PENG S. Effects of mineral dissolution ratios on chemical-dissolution front instability in fluid-saturated porous media [J]. Transport in Porous Media, 2010, 82: 317-335.

[39] HOBBS B E, ORD A. Theoretical analyses of the effects of solute dispersion on chemical-dissolution front instability in fluid-saturated porous media [J]. Transport in Porous Media, 2010, 84: 629-653.

[40] SCHAUBS P M, ZHAO C. Numerical models of gold-deposit formation in the Bendigo-Ballarat Zone, Victoria [J]. Australian Journal of Earth Sciences, 2002, 49: 1077-1096.

[41] SORJONEN-WARD P, ZHANG Y. Numerical modelling of orogenic processes and gold mineralisation in the southeastern part of the Yilgarn Craton, Western Australia [J]. Australian Journal of Earth Sciences, 2002, 49: 935-964.

[42] ZHANG Y, HOBBS B E, ORD A, BARNICOAT A. The influence of faulting on host-rock permeability, fluid flow and ore genesis of gold deposits: A theoretical 2D numerical model [J]. Journal of Geochemical Exploration, 2003, 79: 279-284.

[43] JU M, DAI T, YANG J. Finite element modeling of pore-fluid flow in the Dachang ore district, Guangxi, China: Implications for hydrothermal mineralization [J]. Geoscience Frontiers, 2011, 2: 463-474.

[44] HORNBY P, ORD A, PENG S. Numerical modelling of fluids mixing, heat transfer and non-equilibrium redox chemical reactions in fluid-saturated porous rocks [J]. International Journal for Numerical Methods in Engineering, 2006, 66: 1061-1078.

[45] HORNBY P, PENG S, LIU L. Theoretical and numerical analyses of pore-fluid flow patterns around and within inclined large cracks and faults [J]. Geophysical Journal International, 2006, 166: 970-988.

[46] ZHAO C, HOBBS B E, ORD A. Theoretical analyses of nonaqueous phase liquid dissolution-induced instability in two-dimensional fluid-saturated porous media [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2010, 34: 1767-1796.

[47] ZHAO C, HOBBS B, ALT-EPPING P. Modeling of ore-forming and geoenvironmental systems: Roles of fluid flow and chemical reaction processes [J]. Journal of Geochemical Exploration, 2014, 144: 3-11.

[48] ZHAO C, HOBBS B E, ORD A. A new alternative approach for investigating acidization dissolution front propagation in fluid-saturated rocks [J]. Science China: Technological Sciences, 2017, 60: 1197-1210.

[49] ZHAO C, WANG Y, MUHLHAUS H B, LIN G, ORD A, LIN G. Finite element modelling of reactive fluids mixing and mineralization in pore-fluid saturated hydrothermal/sedimentary basins [J]. Engineering Computations, 2002, 19: 364-387.

[50] POULET T, REGENAUER-LIEB K. Numerical modeling of toxic nonaqueous phase liquid removal from contaminated groundwater systems: Mesh effect and discretization error estimation [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2015, 39: 571-593.

[51] REID L B, REGENAUER-LIEB K, POULET T. A porosity-gradient replacement approach for computational simulation of chemical-dissolution front propagation in fluid-saturated porous media including pore-fluid compressibility [J]. Computational Geosciences, 2012, 16: 735-755.

[52] REID L B, REGENAUER-LIEB K. Some fundamental issues in computational hydrodynamics of mineralization: A review [J]. Journal of Geochemical Exploration, 2012, 112: 21-34.

[53] LIU Y, DAI T, XIA S, TIAN H. Computational simulation of iron ore-forming processes in the Caiyuanzi siderite ore district, Guizhou, China [J]. Journal of Geochemical Exploration, 2015, 158: 155-167.

[54] HOBBS B E, REGENAUER-LIEB K, ORD A. Computational simulation for the morphological evolution of nonaqueous phase liquid dissolution fronts in two-dimensional fluid-saturated porous media [J]. Computational Geosciences, 2011, 15: 167-183.

[55] PENG S, HOBBS B E, ORD A. Particle simulation of spontaneous crack generation associated with the laccolithic type of magma intrusion processes [J]. International Journal for Numerical Methods in Engineering, 2008, 75: 1172-1193.

[56] PENG S, LIU L, HOBBS B E, ORD A. Computational simulation of convective flow in the Earth crust under consideration of dynamic crust-mantle interactions [J]. Journal of Central South University of Technology, 2011, 18: 2080-2084.

[57] ZHAO C, HOBBS B E, ORD A. Theoretical analyses of chemical dissolution-front instability in fluid-saturated porous media under non-isothermal conditions [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2015, 39: 799-820.

[58] ZHAO C, HOBBS B E, ORD A. Computational simulation of chemical dissolution-front instability in fluid-saturated porous media under non-isothermal conditions [J]. International Journal for Numerical Methods in Engineering, 2015, 102: 135-156.

[59] ZHAO C. Physical and chemical dissolution front instability in porous media: Theoretical analyses and computational simulations [M]. Berlin: Springer, 2014.

[60] POULET T, REGENAUER-LIEB K. Replacement of annular domain with trapezoidal domain in computational modeling of nonaqueous-phase-liquid dissolution-front propagation problems [J]. Journal of Central South University, 2015, 22: 1841-1846.

[61] WANG S, FENG C Y, BO H X, JIANG J H. Characteristics and genesis of mineral compostion of Kaerqueka skarn type copper polymetallic deposit in Qimantage Area, Qinghai Province [J]. Acta Mineralogica Sinica, 2009, 29: 483-484.

[62] LI D S, ZHANG W Q, TIAN C S, YAN C, WANG L J, JING X. Discussion on the metallogenic characteristics and ore-prospecting methods of Qimantage region, Qinghai Province [J]. Northwestern Geology, 2013, 46: 131-141. (in Chinese)

[63] XUE P. Geological characteristics and ore-finding prospect analysis of Yelasai Copper Mine in Kaerqueka [J]. West-China Exploration Engineering, 2007, 19: 68-71. (in Chinese)

[64] LI D X, FENG C Y, ZHAO Y M, LI Z F, LIU, J N, XIAO Y. Mineralization and alteration types and skarn mineralogy of Kaerqueka copper polymetallic deposit in Qinghai Province [J]. Journal of Jilin University, 2011, 41: 1818-1830. (in Chinese)

[65] LIANG H, ZHANG S, ZHANG H, CHEN Z, ZHANG Z, CHANG, Y. Geological characteristics and metallogenic model of Kaerqueka copper polymetallic deposit in Geermu of Qinghai [J]. Mineral Resources and Geology, 2015, 29: 7-13. (in Chinese)

[66] SONG Z B, JIA Q Z, ZHANG Z Y, HE S Y, CHEN X Y, QUAN S C, LI Y Z, ZHANG Y L, ZHANG X F. Study on geological feature and origin of Yemaquan Fe-Cu deposit in Qimantage Area, Eastern Kunlun [J]. Northwestern Geology, 2010, 43: 209-217. (in Chinese)

[67] CUNDALL P A. A microcomputer program for modeling large-strain plasticity problems [J]. Numerical Methods in Geomechanics (Innsbruck 1988), 1988, 1: 2101-2108.

[68] LIU L, ZHAO Y. Coupled geodynamics in the formation of Cu skarn deposits in the Tongling-Anqing district, China: Computational modeling and implications for exploration [J]. Journal of Geochemical Exploration, 2010, 106: 146-155.

[69] POULET T, REGENAUER-LIEB K, HOBBS B E. Computational modeling of moving interfaces between fluid and porous medium domains [J]. Computational Geosciences, 2013, 17: 151-166.

[70] ALT-EPPING P, ZHAO C. Reactive mass transport modelling of a three-dimensional vertical fault zone with a finger-like convective flow regime [J]. Journal of Geochemical Exploration, 2010, 106: 8-23.

[71] MUMM A S, BRUGGER J, ZHAO C, SCHACHT U. Fluids in geological processes—The present state and future outlook [J]. Journal of Geochemical Exploration, 2010, 106: 1-7.

[72] WALSHE J L, MUHLHAUS H B, ORD A. Finite element modeling of fluid-rock interaction problems in pore-fluid saturated hydrothermal/sedimentary basins [J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190: 2277-2293.

[73] ORD A, HOBBS B E, ZHAO B C. Analytical solutions of nonaqueous-phase-liquid dissolution problems associated with radial flow in fluid-saturated porous media [J]. Journal of Hydrology, 2013, 494: 96-106.

[74] ORD A, HOBBS B E, ZHAO B C. Effects of medium permeability anisotropy on chemical-dissolution front instability in fluid-saturated porous media [J]. Transport in Porous Media, 2013, 99: 119-143.

[75] ZHAO C. Dynamic and transient infinite elements: Theory and geophysical, geotechnical and geoenvironmental applications [M]. Berlin: Springer, 2009.

[76] ORD A, HOBBS B E, ZHAO B C. Effects of domain shapes on the morphological evolution of nonaqueous-phase-liquid dissolution fronts in fluid-saturated porous media [J]. Journal of Contaminant Hydrology, 2012, 138: 123-140.

[77] CHENG Yong-sheng. Geochemistry of intrusive rock in Dachang tin-polymetallic ore field, Guangxi, China: Implications for petrogenesis and geodynamics [J]. Transactions of Nonferrous Metals Society of China, 2015, 25: 284-292.

[78] HUA Xiao-ming, ZHENG Yong-fei, XU Qian, LU Xiong-gang, CHENG Hong-wei, ZOU Xing-li, SONG Qiu-shi, NING Zhi-qiang. Interfacial reactions of chalcopyrite in ammonia-ammonium chloride solution [J]. Transactions of Nonferrous Metals Society of China, 2018, 28: 556-566.

[79] XU T, SONNENTHAL E, SPYCHER N, PRUESS K. TOUGHREACT user’s guide: A simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media (V1.2.1) [R]. Berkeley, CA (US): Ernest Orlando Lawrence Berkeley National Laboratory, 2008.

青海卡尔却卡多金属矿床热液与化学反应耦合成矿过程的数值模拟

邹艳红^1,2，刘 尧^1,3，潘 勇⁴，阳宽达^1,2，戴塔根^1,2，毛先成^1,2，赖建清^1,2，田海龙⁵

1. 中南大学 计算地球科学研究中心，长沙 410083；

2. 中南大学 地球科学与信息物理学院 有色金属成矿预测与地质环境监测教育部重点实验室，长沙 410083；

3. 湖南大学 嵌入式系统及网络实验室，长沙 410082；

4. 湖南文理学院 资源环境与旅游学院，常德 415000；

5. 吉林大学 地下水资源与环境教育部重点实验室，长春 130021

摘  要：青海卡尔却卡多金属矿床位于祁漫塔格成矿带，是一个典型的矽卡岩多金属矿床，该矿床形成的条件及详细动力学过程一直是地球科学研究工作中的重要课题。利用有限差分法对该地区的黄铜矿成矿过程进行数值建模，并在数值模拟中近似地考虑野外地质特征、成矿及其相关的地球化学条件。着重采用现代成矿理论定量地考虑卡尔却卡多金属矿床中黄铜矿的成矿化学反应过程。相关数值结果表明：含矿热液流动是该地区矿化的关键控制因素，而温度梯度是孔隙热液流体流动的主要驱动力。在卡尔却卡多金属矿床中，黄铜矿成矿的温度为250~350 °C。这些相关的计算结果已通过地质勘查成果得到验证。由此表明：采用新兴计算地球科学中的成矿化学反应与有限差分模型耦合模拟方法有利于提高对该地区成矿过程的认识。

关键词：数值模拟；热传递；化学反应；热液成矿；卡尔却卡多金属矿床

(Edited by Wei-ping CHEN)

Foundation item: Project (2017YFC0601503) supported by the National Key R&D Program of China; Projects (41872249, 41472302, 41772348) supported by the National Natural Science Foundation of China

Corresponding author: Yao LIU; Tel: +86-18608405455; E-mail: liuyao216@163.com

DOI: 10.1016/S1003-6326(18)64925-8