J. Cent. South Univ. (2019) 26: 501-514

DOI: https://doi.org/10.1007/s11771-019-4022-x

Finite element modeling of convective pore-fluid flow in fluid-saturated porous rocks within upper crust: An overview

ZHAO Chong-bin(赵崇斌)¹, Bruce HOBBS², Alison ORD²

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

2. School of Earth and Environment, the University of Western Australia, Crawley, WA 6009, Australia

Central South University Press and Springer-Verlag GmbH Germany, part of Springer Nature 2019

Abstract: Convective pore-fluid flow (CPFF) plays a critical role in generating mineral deposits and oil reservoirs within the deep Earth. Therefore, theoretical understanding and numerical modeling of the thermodynamic process that triggers and controls the CPFF are extremely important for the exploration of new mineral deposits and underground oil resources. From the viewpoint of science, the CPFF within the upper crust can be treated as a kind of thermodynamic instability problem of pore-fluid in fluid-saturated porous media. The key issue of dealing with this kind of problem is to assess whether a nonlinear thermodynamic system under consideration is supercritical. To overcome limitations of using theoretical analysis and experimental methods in dealing with the CPFF problems within the upper crust, finite element modeling has been broadly employed for solving this kind of problem over the past two decades. The main purpose of this paper is to overview recent developments and applications of finite element modeling associated with solving the CPFF problems in large length-scale geological systems of complicated geometries and complex material distributions. In particular, two kinds of commonly-used finite element modeling approaches, namely the steady-state and transient-state approaches, and their advantages/disadvantages are thoroughly presented and discussed.

Key words: convective flow; steady-state approach; transient-state approach; numerical modeling; upper crust; porous rock

Cite this article as: ZHAO Chong-bin, Bruce HOBBS, Alison ORD. Finite element modeling of convective pore-fluid flow in fluid-saturated porous rocks within upper crust: An overview [J]. Journal of Central South University, 2019, 26(3): 501–514. DOI:https://doi.org/10.1007/s11771-019-4022-x.

1 Introduction

Convective pore-fluid flow (CPFF) may occur in a wide range of geological environments, so that it can cause the occurrence of many geological phenomena, for example, the formation of hot spots (including hot springs), mineral deposits and oil reservoirs within the upper crust of the Earth (UCE). From the viewpoint of science, the CPFF is a most widespread and efficient agent that plays a major role in the evolution of geological environments by transporting mass and thermal energy. It can significantly change the geological environments suddenly because of its emerging behavior (in the complex Earth system) or more slowly change over time before it reaches a steady state. Therefore, theoretical understanding and numerical modeling of the thermodynamic process that triggers and controls convective pore-fluid flow are extremely important for exploring new mineral deposits and oil reservoirs, as well as for exploiting of hot spring water within the UCE [1–7]. Due to ever-increasing demand on mineral and energy resources, the mankind is confronting a sharp supply-demand contradictory because the existing mineral and energy resources are exhausted in the foreseeable future. With a powerful drive to understand ore- forming and mineralization mechanisms, extensive theoretical and numerical studies on the controlling thermodynamic process of CPFF within the UCE have been already carried out over the past two decades or so [8, 9].

The specific scientific problem associated with the CPFF within the UCE belongs essentially to a kind of thermodynamic instability problem of pore-fluid in the fluid-saturated porous medium (FSPM) [8, 9]. From the viewpoint of continuum mechanics, we can treat the crustal rock as a porous medium. The pore-fluid (such as water), which resides in the pore space of the porous medium, is initially in a rest state, so that no pore-fluid flow initially occurs within the UCE. When the flat bottom of the upper crust is heated due to the occurrence of a thermal event in the geological history, the initial temperature gradient is changed, so that a new gradient of temperature is formed within the UCE. Since the temperature in the lower part of the upper crust is higher than that in the upper part, there is an upward buoyant force exerted on the pore-fluid. However, if the pore-fluid intends to move upwards, it will be exposed to a resistant force resulting from both its viscosity and the roughness of the flow channels, which are usually generated by connected pores in the porous medium. Because the upward buoyant force is proportional to the pore-fluid temperature, there exists a critical state, in which the upward buoyant force exerted on the pore-fluid is straightforwardly equal to the viscous resistant force of the pore-fluid. This indicates that the thermodynamic system associated with the CPFF can be in one of the following states: a supercritical state, a critical state and a subcritical state. From the viewpoint of complex system science, the loss of stability is called the emerging behavior of a thermodynamic system, while the critical state is called the emerging point of the thermodynamic system. Since this kind of thermodynamic instability of pore-fluid in the FSPM was first investigated by LAPWOOD [10] as well as HORTON et al [11], it was previously named as the Horton-Rogers-Lapwood instability problem [8]. Nevertheless, this kind of thermodynamic instability problem is the counterpart of the Rayleigh-Benard instability problem [12], which is another well-known kind of thermodynamic instability problem involved in the study of mantle convection. The similarity between the Horton-Rogers-Lapwood instability problem and the Rayleigh-Benard instability problem is that the upward buoyant force is the driving force for both the problems, while the difference between them is that the Horton-Rogers-Lapwood instability problem deals with the CPFF within the UCE consisting of porous rocks, but the Rayleigh-Benard instability problem deals with convective viscous- fluid flow within the mantle that is commonly treated as viscous fluid of high viscosity, compared to water. Due to this similarity, the critical Rayleigh number of two different expressions is widely used to assess whether or not the thermodynamic system under consideration is unstable.

By following the original work of LAPWOOD [10] as well as HORTON et al [11], extensive theoretical and experimental studies have been conducted to investigate the CPFF in the FSPM over the past several decades [9, 13–24]. Consequently, the following main conclusions have already been made in the recent years: 1) if a thermodynamic system has either a supercritical or a critical Rayleigh number, then the CPFF can take place within the upper crust consisting of fluid- saturated porous rocks, in which hydrostatic pressure gradient is predominant; 2) If lithostatic pressure gradient is predominant in fluid-saturated porous rocks, then the CPFF cannot take place in the porous rocks unless the bottom boundary of a thermodynamic system is permeable to both pore-fluid and heat flux; 3) Either permeability heterogeneity or thermal conductivity heterogeneity of porous rocks can have significant effects on the CPFF patterns in a supercritical thermodynamic system; 4) CPFF can also occur in three- dimensional large fault zones within the UCE, if the Rayleigh number of a fault system is supercritical.

However, theoretical analysis and experimental methods have some limitations in dealing with CPFF within the UCE that is comprised of porous rocks. For example, it is often very difficult to use them to consider realistic geological systems of complicated geometrical shapes and complex rock (material) distributions within the UCE. With the rapid development of both modern computer technology and computational geosciences, numerical modeling has become very popular in solving many geophysical

and geological problems, which cannot be solved by either the traditionally-used theoretical method or the existing experimental method because of the involvement of both the large timescale and length-scale. Due to the complex and complicated nature of geoscience problems, numerical modeling has found wider and wider applications in almost all fields of geosciences [25–35]. This has led to a significant change in the research methodology of traditional geosciences [36, 37]. Consequently, previously-unsolvable geophysical and geological problems can be nowadays solved through using finite element modeling and other numerical simulation methods. Therefore, it is necessary to discuss the recent development and application of finite element modeling in dealing with large-scale CPFF problems in realistic geophysical and geological systems of complicated geometrical shapes and complex rock (material) distributions within fluid-saturated porous rocks.

2 Brief mathematical description of CPFF problems within upper crust

The CPFF problem within the UCE can be commonly treated as a coupled nonlinear problem involving heat transfer and pore-fluid flow processes in fluid-saturated (rigid) porous rocks. When the axis directions of a rectangular coordinate system are appropriately chosen, the crustal rock may be viewed as an orthotropic porous medium. For example, one axis of the rectangular coordinate system can be parallel to the vertical direction, while another axis is aligned in the horizontal direction for a two-dimensional problem or other two axes are aligned in the horizontal plane for a three-dimensional problem. For the sake of conciseness, only the related mathematical formulations of a two-dimensional problem are briefly given in this overview. Based on the mass conservation of the pore-fluid, together with the appropriate consideration of corresponding scientific principles, the governing mathematical equations of the considered problem can be expressed in the two-dimensional case as follows [9]:

                (1)

                      (2)

                     (3)

     (4)

                     (5)

,

                        (6)

where f is the crustal rock porosity; p is the pore-fluid pressure in the pores of the crustal rock; u and v are the pore-fluid velocity components in the horizontal and vertical directions; ρ_f is the pore-fluid density; T is the crustal rock temperature; ρ_f0 is the pore-fluid reference density; T₀ is the crustal rock reference temperature; μ is the pore-fluid dynamic viscosity; K_y and K_x are the permeability of the crustal rock in the vertical and horizontal directions; Q_f is the pore-fluid source term (representing the addition/generation of the pore fluid per unit volume of the crustal rock); ρ_s is the solid matrix density of the crustal rock; and are the pore-fluid thermal conductivity coefficients in the vertical and horizontal directions, respectively; c_pf and c_ps are the pore-fluid specific heat and solid matrix specific heat in the crustal rock; β_T is the pore-fluid thermal volume expansion coefficient; Q_T is the crustal rock heat source (representing the addition/generation of heat per unit volume of the crustal rock); and are the crustal (dry) rock thermal conductivity coefficients in the vertical and horizontal directions, respectively; g_y and g_x are the components of gravity acceleration in the y and x directions, respectively. Generally, it is assumed that gravity is opposite to the positive direction of the y axis.

Note that Eq. (1) represents the pore-fluid continuous equation, which is used for describing the pore-fluid conservation; Eqs. (2) and (3) represent Darcy’s law, which are used, respectively, for describing the pore-fluid momentum conservation in the horizontal and vertical directions; Eq. (4) represents the energy equation, which is used for describing the energy conservation within the crustal rock. Because the pore-fluid density depends on temperature (see Eq. (5)), the temperature and pore-fluid flow fields are strongly coupled through the energy equation and pore-fluid velocity (see Eq. (1) and Eq. (4)).

Specifically, the pore-fluid (i.e. water) dynamic viscosity can be written into a function of temperature below [38]:

     (7)

                            (8)

where T₀ is the pore-fluid reference temperature in the crustal rock; μ₀ is the pore-fluid reference dynamic viscosity in correspondence to the reference temperature, T₀.

The pore-fluid density may depend on salinity, temperature and pressure. In some cases, pore-fluid such as water may be less compressible, so that the pore-fluid density becomes almost independent of pore-fluid pressure. If the pore-fluid density variation is small or moderate, then the commonly-used Boussinesq approximation [8, 9] is valid. This means that in such a case, the pore-fluid density is also independent of salinity. Based on the Boussinesq approximation, which stated that the pore-fluid (such as water) density can be treated as a constant in both the pore-fluid continuity equation and the energy equation, it is possible to simplify Eqs. (1) and (4) as follows:

                           (9)

    (10)

To justify the validity of using the Boussinesq approximation in theoretical analysis and finite element modeling of CPFF problems within the UCE, we need to examine the difference between the original form (i.e. Eq. (1)) and the simplified form (i.e. Eq. (9)) of the pore-fluid continuity equation. If we insert Eq. (5) into Eq. (1), then we can obtain the following equation:

 (11)

This equation clearly shows that the original form (i.e. Eq. (1)) differs from the simplified form (i.e. Eq. (9)) of the pore-fluid continuity equation in the following way:

     (12)

Equation (12) is called the term of the extended Boussinesq approximation. This means that an additional term should be considered for incorporating temperature-dependent effects in the simplified pore-fluid continuity equation. Note that if the considered hydrothermal system reaches a steady state, then the temporal changes can vanish in temperature, so that the first term of Eq. (12) becomes unimportant. Similarly, if the pore-fluid flow is orthogonal to the temperature gradient, then the second term of Eq. (12) can also vanish. On the other hand, because of relatively small values of the pore-fluid thermal volume expansion coefficient, which is typically equal to around 2.07×10^{–4 °C–1} for water, the overall value of the term of the extended Boussinesq approximation is relatively small in the original pore-fluid continuity equation. For this reason, the Boussinesq approximation is widely accepted in theoretical analysis and finite element modeling of CPFF problems within the CUE [8, 9, 25].

If both the source term of the pore-fluid and the heat source term in the crustal rock are neglected, then the simplified pore-fluid continuity equation and energy equation can be further expressed below:

                             (13)

                    (14)

Finally, if we insert Eq. (5) into Eqs. (2) and (3), then we can obtain the following equations:

    (15)

    (16)

When a hydrothermal system reaches a steady state with T/t=0, Eqs. (13)–(16) form the fundamental governing mathematical equations of the Horton-Rogers-Lapwood instability problem associated with the CPFF within the UCE [8, 9]. As demonstrated from a previous study [39], if an upper crust has a thickness of 10 km, then it takes around 7.79 million years, so that the whole upper crust could be warmed up from the bottom when the pure heat conduction process is only considered. This means that if a thermal event can last long enough (e.g. longer than 7.79 million years) in the geological history, then it is possible to enable the hydrothermal system to attain a steady state within the UCE. This is the main reason why the Horton-Rogers-Lapwood instability problem is commonly investigated in theoretical analysis and finite element modeling of the CPFF within the UCE [8, 9].

3 Finite element modeling of CPFF problems within upper crust

From the viewpoint of physics, CPFF occurrence within fluid-saturated porous rocks is an emerging behaviour of the complex Earth system. If a hydrothermal system has a Rayleigh number that reaches or exceeds the corresponding critical Rayleigh number, then the CPFF can occur instantaneously. This emerging phenomenon is companied with a dissipative structure, namely the CPFF pattern in the thermodynamic instability system. More importantly, such a dissipative structure has a steady-state characteristic, so that it is time-independent if we can maintain the thermodynamic properties to be unchanged in the thermodynamic system. Therefore, CPFF can be simulated computationally and analyzed mathematically through the steady-state analysis of a hydrothermal system. On the basis of this fact, the steady-state approach was commonly employed to investigate theoretically the convective pore-fluid flow instability in hydrothermal systems within fluid-saturated porous rocks [8, 9]. To use the steady-state approach efficiently, the progressive asymptotic approach procedure (PAAP) [27] is developed to numerically simulate the steady-state CPFF within the UCE. Compared with the transient-state approach, the primary advantages of using the steady-state approach are as follows:1) since the CPFF problem is mathematically treated as a boundary-value problem, the temperature-field build-up process (that is a transient process) within the UCE is skipped, so that the computational simulation time can be greatly reduced in the finite element modeling of CPFF within the UCE; 2) since the initial condition, which is usually associated with the past in the geological history, is not involved in a steady-state problem, some uncertainty resulting from using the initial condition can be removed from the finite element modeling of CPFF. The primary disadvantage in using the steady-state approach is that the temperature-field build-up process within the UCE cannot be considered in finite element modeling of CPFF.

However, if one is interested in simulating the temperature-field build-up process within the UCE, then the transient-state approach needs to be used to computationally simulate CPFF problems within the UCE. Compared with the steady-state approach, the primary advantage in using the transient-state approach is that since CPFF problem mathematically belongs to an initial-value problem, we can simulate both temperature-field build-up and pore-fluid flow-field build-up processes within the UCE. The primary disadvantage in using the transient-state approach is that it is necessary to use more computer efforts because many time steps need to be used in finite element modeling of CPFF within the UCE, unless the problem itself eventually belongs to a transient one.

3.1 Steady-state approach based on finite element method

To briefly illustrate how the finite element method (FEM) [40] can be used to deal with large-scale CPFF problems within fluid-saturated heterogeneous rocks, a two-dimensional problem is considered to derive the related mathematical formulas. For this purpose, the distribution of the four unknown variables (i.e. pore-fluid pressure, temperature, vertical and horizontal components of the velocity) associated with a large-scale CPFF problem can be described within a finite element in the following manner:

                           (17)

                           (18)

                           (19)

                          (20)

where v^e(x, y) and u^e(x, y) are the distribution fields of vertical and horizontal components of the pore-fluid velocity within the finite element; p^e(x, y) and T ^e(x, y) are the distribution fields of pore-fluid pressure and temperature within the finite element, respectively; v^e and u^e are the nodal column vectors of vertical and horizontal components of the pore-fluid velocity for the finite element; T ^e and p^e are the nodal column vectors of the pore-fluid temperature and pressure for the finite element, respectively; N is the finite element shape function matrix.

According to the conventional procedure of the FEM [9, 40], the discretized governing mathematical equations of the CPFF problem can be expressed, in the finite element sense, as follows:

          (21)

          (22)

          (23)

                        (24)

Note that in a matrix form, we can rewrite Eqs. (21)–(24) as:

     (25)

where A_x^e, A_y^e, B_x^e, B_y^e, C_x^e, C_y^e, E^e and M^e are the property matrices of the finite element; F_x^e, F_y^e and G^e are the nodal load vectors of the finite element. They can be derived and expressed below:

,

              (26)

,

              (27)

,

                       (28)

,

                    (29)

,

                    (30)

               (31)

,

           (32)

,

                         (33)

where S and A are the boundary length and area of the finite element; q_n is the heat flux in the normal direction of the element boundary.

To eliminate the nodal vector of pore-fluid pressure, p^e, we can use the penalty finite element approach (PFEA) [40] to derive the following equation:

                   (34)

where ε is a small parameter used in the penalty finite element approach. The physical meaning of this parameter lies in the fact that we assume that the pore-fluid is slightly compressible in the computational simulation.

It is possible to straightforwardly express Eq. (34) in the following way:

              (35)

If we insert Eq. (35) into Eq. (25), then it is also possible to obtain the following discretized equation at the elemental level:

         (36)

,

,

,

               (37)

Through assembling the discretized equations of all the finite elements appropriately, the global discretized equation of the whole system can be obtained.

   (38)

where Q₁₁, Q₁₂, Q₂₁, Q₂₂, B_x, B_y and E(u, v) are the global property sub-matrices of the whole system; K_global is the global stiffness matrix of the whole system; F_x, F_y and G are the global nodal load sub-vectors of the whole system.

Since E(u, v) depends on the pore-fluid velocity, we need to solve Eq. (38) simultaneously and iteratively. In addition, the PAAP [9] should be used to solve this equation, so as to obtain the steady-state convective solutions within the UCE.

Note that the primary advantage in using the PFEA is that because the pore-fluid pressure (as a variable) is eliminated from the final finite element formulation (i.e. Eq. (38)), the requirement for the memory of a computer can be greatly reduced, especially when the simultaneous solver, in which the final finite element formulation (i.e. Eq. (38)) of the system is solved for temperature and pore-fluid velocities simultaneously, is used to solve the steady-state CPFF problems within the UCE.

The steady-state approach, which is based on the FEM, has been employed for solving many CPFF and related mineralization problems within the UCE. For example, this approach has been successfully used for understanding the pore-fluid convection and related ore-forming mechanisms of Au deposit in Western Australia Yilgarn Craton [5, 36], Queensland Mount Isa Pb–Zn deposit, Australia [3, 41], Cu–Au mineralization in the New Guinea [2], New South Wales Broken Hill Pb–Zn deposit, Australia [42, 43], and other generic ore-forming systems [44].

3.2 Transient-state approach based on finite difference and element methods

In this situation, the FEM is still used for discretizing the computational space. However, the finite difference method (FDM) needs to be used for discretizing the computational time. Again, for the sake of conciseness, a two-dimensional problem is also considered to derive the related mathematical formulation.

If we insert Eqs. (15) and (16) into Eq. (13), then we can obtain the following pore-fluid mass conservation equation:

            (39)

By means of using the conventional procedure of the FEM [9, 40], we can further express Eq. (39) below:

      (40)

where v_n represents the pore-fluid velocity in the normal direction of the element boundary.

Similarly, we can express Eq. (40) in the following manner:

         (41)

where R_x^e, R_y^e, H_x^e and H_y^e are the property matrices of the finite element; J^e is the nodal load vector of the finite element. They can be derived and expressed as follows:

,

                   (42)

,

           (43)

                          (44)

By assembling the discretized equations of all the finite elements adequately, we can express the discretized pore-fluid mass conservation equation of the whole system below:

               (45)

where R_x, R_y, H_x and H_y are the property matrices of the whole system; J is the nodal load vector of the whole system.

Note that Eq. (45) is the global finite element formulation of the pore-fluid mass conservation equation. According to Eq. (25), we can express the global finite element formulation of Darcy’s law (which can be also viewed as the momentum conservation equation) and the energy conservation equation below:

                     (46)

                     (47)

                     (48)

where A_x, A_y, B_x, B_y, C_x, C_y, E and M are the property matrices of the whole system; F_x, F_y and G are the nodal load vectors of the whole system. β can be expressed as:

                   (49)

Based on the segregated algorithm [37], the sequential solver, in which Eqs. (45)–(48) are solved separately, sequentially and iteratively, can be used to solve the transient-state CPFF problems within the UCE.

If the previous computational time-step is t_i, then the current computational time-step is equal to t_i+Δt. Note that Δt represents the computational time-step length hereafter. According to this definition, we can express Eqs. (45)–(48) at the current computational time-step as follows:

          (50)

             (51)

             (52)

      (53)

where p_t+Δt represents the global pore-fluid pressure vector of the whole system at the current computational time-step; T_t and T_t+Δt represent the global temperature vectors of the whole system at the previous computational time-step and current computational time-step respectively; u_t+Δt and v_t+Δt represent the x and y components of the global velocity vectors of the whole system at the current computational time-step.

As mentioned above, we can solve Eqs. (50)–(53) iteratively, separately and sequentially for the global pore-fluid pressure vector, x and y components of the global velocity vectors, as well as global temperature vector at the current time-step. When we solve Eq. (50), we still do not know the global temperature vector at the current time-step. Likewise, when we solve Eqs. (51) and (52), we also do not know the global temperature vector at the current time-step. This clearly means that we need to use an iteration scheme to solve these four equations sequentially, because they are fully coupled in the mathematical sense. At the first iteration step, we can use the global temperature vector, which is already determined at the previous time-step, as a reasonable estimate for the global temperature vector at the current time-step when we solve Eq. (50) for the global pore-fluid pressure vector. In a similar way, we can use the global temperature vector at the previous time-step as a reasonable estimate for the global temperature vector at the current time-step when we solve Eqs. (51) and (52) for the x and y components of the global velocity vectors of the whole system respectively. The resulting approximate global velocity vectors can be used to upgrade the property matrix, E(u, v)_t+Δt, when we solve Eq. (53) for the global temperature vector. At the second iteration step, it is possible to follow the exactly same procedure as we have already adopted at the first iteration step. This allows us to establish a convergence criterion, which can be used to check whether the solution error between two time steps is within the allowable error range, after the second iteration step. Thus, after the second iteration step, it is possible to use this convergence criterion to check whether or not a convergence solution is obtained at the present time-step. If we do not obtain a convergence solution, then we need to repeat the iteration at the present time-step. On the contrary, if we obtain a convergence solution at the present time-step, then we should go to the next time-step until we reach the final time-step.

The transient-state approach, which is essentially based on the FEM and FDM, has also been used to solve many CPFF and related ore-forming problems within the UCE. For instance, JU et al [7] used this approach to simulate the CPFF in the Guangxi Dachang ore field, China. ALT- EPPING et al [35] used it to conduct numerical modeling of a reactive mass transport problem within a vertical fault zone of three-dimensional geometry, in which a finger-like CPFF is the main agent of reactive mass transport. YANG [31] used this approach to carry out the numerical simulation of three-dimensional CPFF in the Northern Australia McArthur basin. In addition, YANG et al [32, 33] also used it to simulate the basin-scale CPFF in the Northern Australia Mount Isa basin. More recently, the transient-state approach has been applied to simulate CPFF problems in horizontal layers [45], fracture networks embedded in a low-permeability matrix [46] and Saskatchewan Athabasca Basin, Canada [47].

4 Application examples

The first application example is designed to investigate the robustness and correctness of the proposed steady-state approach [2, 3, 9]. For this particular purpose, we consider a simplified problem, which contains fluid-saturated homogeneous rocks just in a small portion of the whole UCE. The problem associated with this application example is solved on the basis of directly using dimensional governing mathematical equations of the nonlinear hydrothermal system. The computational domain of this application example has both a length of 110 km and a depth of 13 km, respectively, as shown in Figure 1. It is discretized into 3956 quadrilateral finite elements. Note that the length versus depth ratio of this computational domain, which is greater than 8, can be employed to determine the total number of convective cells in the CPFF system. We can use this ratio to identify whether a CPFF system belongs to a long system, in which multiple CPFF cells should occur. Since the length versus depth ratio of this computational domain is greater than 8, the considered CPFF system belongs to a long system. Theoretically, if the proposed steady-state approach is robust enough, then it should simulate multiple CPFF cells in the finite element model. This forms the fundamental idea behind using an alternative way to verify the proposed steady-state approach. In addition, when the porous rock of the finite element model is assumed to be homogeneous and isotropic, the analytical solution for the CPFF problem considered in this application example is available [8, 9], so that this application example can provide a benchmark solution for examining the robustness of the proposed steady-state approach before it is employed for solving any realistic CPFF problems of complicated geometries and heterogeneous rocks. Through comparing the analytical solutions with computer-simulated ones for this application example, we can answer why the proposed steady-state approach is suitable for dealing with steady-state CPFF problems of complicated geometries and heterogeneous rocks in naturally-formed hydrothermal systems, which will be considered in the second application example later.

Figure 1 Finite element model of CPFF problem in fluid-saturated homogeneous rocks

For the purpose of reflecting geophysical and geological reality, we assume that the considered finite element model has a temperature of 20 °C at its top and a temperature of 150 °C at its bottom. Since the bottom is hotter than the top, the positive direction of thermal flux is from its bottom to its top, resulting in a vertical temperature gradient in the considered computational domain. Because of this vertical temperature gradient, the CPFF may occur in a supercritical hydrothermal system within the UCE. On the vertical boundaries, we need to prevent both mass and thermal energy from losing, so that we should apply both impermeable and insulation horizontal-boundary conditions to them. In addition, on either the top or the bottom boundaries, we need to prevent mass from losing, so that we apply impermeable vertical-boundary conditions to them in the considered finite element model.

We used the following parameters in the forthcoming finite element modeling. We assumed that: the reference density and thermal volume expansion coefficient of the pore-fluid are respectively equal to 1000 kg/m³ and 2.07×10^{–4 °C–1}; the dynamic viscosity of the pore-fluid is equal to 10^–3 N·s/m²; the specific heats are equal to 4185 J/(kg·°C) and 815 J/(kg·°C) for the pore-fluid and porous rock respectively; the porosity of the porous rock is equal to 0.1; the thermal conductivity coefficients are respectively equal to 3.35 W/(m·°C) and 0.6 W/(m·°C) for the porous rock and pore-fluid in either vertical or horizontal direction; the permeability of the porous rock is equal to 10^{–14 m2} in the whole considered computational domain.

Figure 2 shows the temperature and pore-fluid velocity distributions associated with the first application example, which belongs to a CPFF problem in fluid-saturated homogeneous rocks. In particular, Figure 2(a) clearly shows that as expected theoretically, multiple CPFF cells occur in the considered CPFF problem of a long system. Therefore, we can demonstrate that the proposed steady-state approach is suitable for simulating realistic CPFF problems. It can be observed that except for at two vertical boundaries, the distribution pattern of CPFF cells has a clear regularity: both the left-hand-side and right-hand- side neighbors of a clockwise CPFF cell are two anti-clockwise CPFF cells, and vice versa. As shown in Figure 2(b), the temperature distribution is obviously localized. Such temperature localization is mainly caused by nonlinear coupling between the temperature and CPFF fields. This further demonstrates the importance of considering temperature-flow coupling effects in realistic CPFF problems.

As the second application example, the steady-state approach is used to simulate the steady-state CPFF within the North West Shelf basin, Australia [48]. In this case, the part of the UCE to be simulated in the finite element model is comprised of different rocks and faults, which are represented by their different values of permeability. Figure 3 shows the geological model of this basin. This geological model is simulated by the finite elements to form a finite element model. It is worth pointing out that we use different length-scales in the vertical and horizontal directions of this figure, so as to display the detailed geological structures in this geological model. In this figure, layer A is used to display the post-rift, while layers B and C are used to display the pre-rift. In addition, layer D was used to display the syn-rift basin fill and layer E is used to display a source rock, in which the crude oil might be generated. Faults are denoted by F. The basin geometry represents a typical graben with footwall uplift relative to the hangingwall of the major faults. The length and depth of the computational domain to be simulated in the second application example are also equal to 110 and 13 km, respectively. This computational domain is discretized into 3672 quadrilateral elements. To facilitate the forthcoming discussions, the rock of layer C is used as the reference rock. Except for using different permeabilities for different layers and faults, other parameters are the same as those used in the first application example. For the sake of comparison, the permeability of the reference rock in layer C is assumed to be K_C=10^{–14 m2}. In order to reflect the permeability heterogeneity of the considered different rocks, we use K_A/K_C=0.2, K_B/K_C=5, K_D/K_C=0.2, K_E/K_C=1.0 and K_F/K_C=10 in the computational domain. Note that generally, we can use symbol K_i to represent the permeability of item i, where i=A, B, C, D, E, F, as shown in Figure 3.

Figure 2 Simulation results of CPFF problem in fluid-saturated homogeneous rocks:

Figure 4 shows the temperature and pore-fluid velocity distributions associated with the second application example, which belongs to a CPFF problem in fluid-saturated heterogeneous rocks. It is observed that the CPFF predominates within the two faults and the more permeable layer (i.e. layer B). However, due to the relative small thickness of the fault, the CPFF cell is not formed directly within the fault itself, even though it can be well formed within the more permeable layer (i.e. layer B) of 4 km thick. This indicates that the heterogeneity of rock permeability plays an important role in controlling the CPFF pattern within the North West Shelf basin. What can be also clearly noted is that the geometry of the basin (including the fault locations) can have considerable effects on the temperature distributions within the geological model of the North West Shelf basin.

Figure 3 Schematic diagram of geological model representing Australia North West Shelf basin (Unit: km)

Figure 4 Simulation results from geological model of Australia North West Shelf basin:

It needs to be pointed out that the chemical dissolution of a porous rock can cause a change in the porosity and permeability of the porous rock, so that it can form a preferential pore-fluid flow channel in the porous rock [49–52]. As a result, the chemical dissolution of a porous rock may affect CPFF in the fluid-saturated porous rock. This means that investigating the effect of the chemical dissolution of a porous rock on the CPFF in the fluid-saturated porous rock becomes an important new topic in the future research.

5 Conclusions

With the strong drive of finding mineral and oil resources in the deep Earth, mathematical and numerical modeling of CPFF within the UCE has attracted considerable attentions over the past two decades. Consequently, both the transient-state and steady-state approaches have been developed and used to solve CPFF problems within the UCE. The primary advantages in using the steady-state approach are that: 1) since CPFF problems are mathematically treated as boundary-value problems, the temperature-field build-up process (that is a transient process) within the upper crust is skipped, so that the computational simulation time is greatly reduced; 2) since the initial condition, which is usually associated with the past in the geological history, is not involved in a steady-state problem, some uncertainty resulting from using the initial condition can be removed from the finite element modeling of CPFF. The primary disadvantage in using the steady-state approach is that the temperature-field build-up process in the UCE cannot be simulated in the finite element modeling of CPFF.

On the contrary, the primary advantage in using the transient-state approach is that since the CPFF problem mathematically belongs to an initial-value problem, both temperature-field build-up and pore-fluid flow-field build-up processes in the UCE can be realistically simulated. However, the primary disadvantage in using the transient-state approach is that it is necessary to use more computer efforts because many time steps need to be used in finite element modeling of CPFF within the UCE, if the problem itself can be in principle treated as a steady-state one.

Due to the nonlinear coupling feature of the CPFF problem, two kinds of solvers have been commonly used. The first kind of solver is called the simultaneous solver, in which the discretized coupling mathematical equations for the temperature, pore-fluid velocity and pore-fluid pressure are solved simultaneously and iteratively. The second kind of solver is called the sequential solver, in which the discretized coupling mathematical equations for the temperature, pore-fluid velocity and pore-fluid pressure are solved in an iterative, separate and sequential manner. The primary advantage in using the simultaneous solver is that it can obtain a convergent solution more quickly, so that we can reduce the total simulation time significantly. The primary disadvantage in using the simultaneous solver is that since the global “stiffness” matrix of the system needs to be used, the requirement for the computer storage is greatly increased. However, the primary advantage in using the simultaneous solver becomes the primary disadvantage in using the sequential solver, while the primary disadvantage in using the simultaneous solver becomes the primary advantage in using the sequential solver.

References

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

[2] GOW P, UPTON P, ZHAO C, HILL K. Copper-gold mineralization in the New Guinea: Numerical modeling of collision, fluid flow and intrusion-related hydrothermal systems [J]. Australian Journal of Earth Sciences, 2002, 49(4): 753–771.

[3] ORD A, HOBBS B E, ZHANG Y, BROADBENT G C, BROWN M, WILLETTS G, SORJONEN-WARD P, WALSHE J, ZHAO C. Geodynamic modelling of the Century deposit, Mt Isa Province, Queensland [J]. Australian Journal of Earth Sciences, 2002, 49(6): 1011–1039.

[4] SCHAUBS P, ZHAO C. Numerical modelling of gold-deposit formation in the Bendigo-Ballarat zone, Victoria [J]. Australian Journal of Earth Sciences, 2002, 49(6): 1077–1096.

[5] SORJONEN-WARD P, ZHANG Y, ZHAO C. Numerical modelling of orogenic processes and mineralization in the south eastern part of the Yilgarn Craton, Western Australia [J]. Australian Journal of Earth Sciences, 2002, 49(6): 935–964.

[6] HARCOUET-MENOU V, GUILLOU-FROTTIER L, BONNEVILLE A, ADLER P M, MOURZENKO V. Hydrothermal convection in and around mineralized fault zones: Insights from two- and three-dimensional numerical modeling applied to the Ashanti belt, Ghana [J]. Geofluids, 2009, 9(2): 116–137.

[7] 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(3): 463–474.

[8] NIELD D A, BEJAN A. Convection in porous media [M]. New York: Springer-Verlag, 1992: 356.

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

[10] LAPWOOD E R. Convection of a fluid in a porous medium [J]. Proceedings of the Cambridge Philosophical Society, 1948, 44(4): 508–521.

[11] HORTON C W, ROGERS F T. Convection currents in a porous medium [J]. Journal of Applied Physics, 1945, 16(6): 367–370.

[12] TURCOTTE D L, SCHUBERT G. Geodynamics: Applications of continuum physics to geological problems [M]. New York: John Wiley & Sons, 1982.

[13] GASSER R D, KAZIMI M S. Onset of convection in a porous medium with internal heat generation [J]. ASME Journal of Heat Transfer, 1976, 98(1): 49–54.

[14] HORNE R N, CALTAGIRONE J P. On the evaluation of thermal disturbances during natural convection in a porous medium [J]. Journal of Fluid Mechanics, 1980, 100(2): 385–395.

[15] BAU H H, TORRANCE K E. Low Rayleigh number thermal convection in a vertical cylinder filled with porous materials and heated from below [J]. ASME Journal of Heat Transfer, 1982, 104(1): 166–172.

[16] KAVIANY M. Thermal convective instabilities in a porous medium [J]. ASME Journal of Heat Transfer, 1984, 106(1): 137–142.

[17] CALTAGIRONE J P, BORIES S. Solutions and stability criteria of natural convective flow in an inclined porous layer [J]. Journal of Fluid Mechanics, 1985, 155: 267–287.

[18] LEBON G, CLOOT A. A thermodynamical modeling of fluid flows through porous media: Application to natural convection [J]. International Journal of Heat and Mass Transfer, 1986, 29(3): 381–390.

[19] PILLATSIS G, TASLIM M E, NARUSAWA U. Thermal instability of a fluid-saturated porous medium bounded by thin fluid layers [J]. ASME Journal of Heat Transfer, 1987, 109(3): 677–682.

[20] BJORLYKKE K, MO A, PALM E. Modelling of thermal convection in sedimentary basins and its relevance to diagenetic reactions [J]. Marine and Petroleum Geology, 1988, 5(4): 338–351.

[21] ALAVYOON F. On natural convection in vertical porous enclosures due to prescribed fluxes of heat and mass at the vertical boundaries [J]. International Journal of Heat and Mass Transfer, 1993, 36(10): 2479–2498.

[22] CHEVALIER S, BERNARD D, JOLY N. Natural convection in a porous layer bounded by impervious domains: From numerical approaches to experimental realization [J]. International Journal of Heat and Mass Transfer, 1999, 42(4): 581–597.

[23] TOURNIER C, GENTHON P, RABINOWICZ M. The onset of natural convection in vertical fault planes: Consequences for the thermal regime in crystalline basements and for heat recovery experiments [J]. Geophysical Journal International, 2000, 140(3): 500–508.

[24] 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(10): 787–799.

[25] PHILLIPS O M. Flow and reactions in permeable rocks [M]. Cambridge: Cambridge University Press, 1991: 286.

[26] RAFFENSPERGER J P, GARVEN G. The formation of unconformity-type uranium ore deposits: Coupled hydrochemical modeling [J]. American Journal of Science, 1995, 295(6): 639–696.

[27] ZHAO C, MUHLHAUS H B, HOBBS B E. Finite element analysis of steady-state natural convection problems in fluid-saturated porous media heated from below [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1997, 21(12): 863–881.

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

[29] 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(7): 1061–1078.

[30] SCHAFER D, SCHAFER W, KINZELBACH W. Simulation of reactive processes related to biodegradation in aquifers: 1. Structure of the three-dimensional reactive transport model [J]. Journal of Contaminant Hydrology, 1998, 31(1, 2): 167–186.

[31] YANG J W. Full 3D numerical simulation of hydrothermal fluid flow in faulted sedimentary basins: Example of the Mcarthur basin, northern Australia [J]. Journal of Geochemical Exploration, 2006, 89(1–3): 440–444.

[32] YANG J W, FENG Z, LUO X, CHEN Y. Three-dimensional numerical modeling of salinity variations in driving basin-scale ore-forming fluid flow: Example from Mount Isa Basin, northern Australia [J]. Journal of Geochemical Exploration, 2010, 106(1–3): 236–243.

[33] YANG J W, LARGE R, BULL S, SCOTT D. 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 [J]. Economic Geology, 2006, 101(6): 1275–1292.

[34] PALUSZNY A, MATTHAI S K, HOHMEYER M. Hybrid finite element–finite volume discretization of complex geologic structures and a new simulation workflow demonstrated on fractured rocks [J]. Geofluids, 2007, 7(2): 186–208.

[35] ALT-EPPING P, ZHAO C. Reactive mass transport modeling of a three-dimensional vertical fault zone with a finger-like convective flow regime [J]. Journal of Geochemical Exploration, 2010, 106(1–3): 8–23.

[36] ZHAO C, 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(7): 947–954.

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

[38] DIERSCH H J G. FEFLOW Reference manual [M]. Berlin: Wasy GmbH, 2002: 225.

[39] ZHAO C, 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.

[40] ZIENKIEWICZ O C. The finite element method [M]. London: McGraw-Hill, 1977: 586.

[41] KUHN M, DOBERTB F, GESSNER K. Numerical investigation of the effect of heterogeneous permeability distributions on free convection in the hydrothermal system at Mount Isa, Australia [J]. Earth and Planetary Science Letters, 2006, 244(3, 4): 655–671.

[42] ZHAO C, HOBBS B E, 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.

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

[44] ZHAO C, SCHAUBS P, HOBBS B E. 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.

[45] MALKOVSKY V I, PEK A A. Onset of fault-bounded free thermal convection in a fluid-saturated horizontal permeable porous layer [J]. Transport in Porous Media, 2015, 110(1): 25–39.

[46] VUJEVIC K, GRAF T. Combined inter- and intra-fracture free convection in fracture networks embedded in a low-permeability matrix [J]. Advances in Water Resources, 2015, 84: 52–63.

[47] PEK A A, MALKOVSKY V I. Linked thermal convection of the basement and basinal fluids in formation of the unconformity-related uranium deposits in the Athabasca Basin, Saskatchewan, Canada [J]. Geofluids, 2016, 16(5): 925–940.

[48] ZHAO C. 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(5): 783–795.

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

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

[51] 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, 2019, 62(1): 163–174.

[52] ZHAO C, HOBBS B, ORD A. Effects of different numerical algorithms on simulation of chemical dissolution-front instability in fluid-saturated porous rocks [J]. Journal of Central South University, 2018, 25(8): 1966–1975.

(Edited by FANG Jing-hua)

中文导读

上地壳内饱水孔隙岩石中孔隙流体对流的有限元模拟综述

摘要：孔隙流体对流在生成地下深部矿产资源和油田过程中起着关键作用。因此，为了探测新的地下深部矿产资源和油田，非常有必要对驱动和控制饱水孔隙岩石中孔隙流体对流的热动力过程进行理论分析和数值模拟。根据科学的观点，上地壳内孔隙流体对流问题可被归结为一类发生在饱水孔隙介质中的热动力非稳定性问题。处理这类科学问题的关键点在于如何评价所考虑的热动力系统是否处于超临界状态。为了克服采用理论分析和实验方法在求解上地壳内孔隙流体对流问题时的局限性，有限元模拟方法已在过去二十多年中发展成为广泛使用的一种有效方法。本文的主要目的是对采用有限元模拟方法求解上地壳内孔隙流体对流问题的发展过程及应用进行综述。所考虑的应用主要涉及采用有限元模拟方法求解具有复杂几何构形和介质材料分布的大尺度地质系统中孔隙流体对流问题。尤其重要的是，本文详细地介绍了两种常用的有限元数值模拟方法，即稳态方法和瞬态方法，并对它们在模拟上地壳内孔隙流体对流问题时的优缺点进行了讨论。

关键词：孔隙流体对流；稳态方法；瞬态方法；数值模拟；上地壳；孔隙岩石

Foundation item: Project(11272359) supported by the National Natural Science Foundation of China

Received date: 2018-04-30; Accepted date: 2018-09-30

Corresponding author: ZHAO Chong-bin, PhD, Professor; Tel: +86-731-88830055; E-mail: chongbin.zhao@iinet.net.au; ORCID: 0000- 0001-7093-2282