Migration of the Industrial Wastewater in Fractured Rock Masses Based on the Thermal-Hydraulic-Mechanical Coupled Model

Industrial wastewater may have a long-time e ﬀ ect on the environment and human life as it goes underground and causes serious pollution continuously. To have a well understanding of the migration of such wastewater is a basic task for industrial wastewater treatment as well as industrial design. To study the migration mechanism of industrial wastewater in rock formations, the governing equations such as mechanics, seepage, heat, and mass transfer are reviewed, referenced, and proposed. The thermal (T)-hydraulic (H)-mechanical (M) coupled model of the multimedia of matrix-fault and matrix-fracture-fault is established. The in ﬂ uence of the fault and the fractures on the pressure distribution and contaminant migration is analyzed. The in ﬂ uence of fault length, width, dip angle, permeability, and temperature of wastewater on contaminant migration is parametrically studied. The following results can be obtained. (1) The fracture quantitively a ﬀ ects the concentration distribution, while the fault dominates the concentration distribution and contaminant migration. (2) The migration of the contaminants can be geometrically divided into 3 zones along the direction of the fault: the saturation zone, the rapid di ﬀ usion zone, and the concentration decrease zone. (3) There is a peak of the concentration along the bottom of the model. The position of the peak is the projection of the endpoint of the fault. (4) The fault length has the most signi ﬁ cant e ﬀ ect on contaminant accumulation. The temperature of the wastewater has the minimum e ﬀ ect on the contaminant accumulation. (5) The accumulation of concentrations can be divided into 2 stages, the slow growth stage (before 20 years) and the rapid growth stage (after 20 years). The main channel of contaminant migration in the slow growth stage is a fault. During the rapid growth stage, the contaminants penetrate through the rock matrix as well as the fault.


Introduction
The great development of the society and economy of China leads the rapid growth in the industry. However, the development of the industry inevitably results in environmental problems as the wastewater discharge cannot be avoided. The wastewater that contains the contaminants will cause long-time pollution and risks for the social life and ecology because it migrates continuously in the underground, such as soil and rock. Such polluted soil or rock thus becomes poi-sonous and harmful and has adverse effects on the environment and human life [1][2][3]. Therefore, the study of the transport process of contaminants in industrial wastewater in rock formations has been a hot topic for scholars in the last decades.
The simultaneous existence of different scales of physical structures in rock formations, such as matrix, fracture networks, and faults, leads to strongly nonhomogeneous behaviors [4][5][6][7]. The contaminants in the migration process thus penetrate from even microscale to macroscale. The flow status also changes from linear Darcy seepage in the rock and natural microfracture network to the complicated nonlinearity in faults. Meanwhile, the transport of contaminants in industrial wastewater in rock formation is the result of a combination of multiple physical fields, including thermal field, hydraulic field, and stress field. The fluid pressure changes in the hydraulic field due to the fluid flow behavior and affects the stress distribution in the rock formation, i.e., it influences the stress field. Also, the thermal characteristics are influenced by the change of velocity of the flow. The change of porosity and permeability of the rock formation due to the change of the stress field also affects the hydraulic field. The change of temperature in the thermal field influences the hydrodynamic viscosity coefficient and the fluid density, which affects the hydraulic field. Thermal stress has a consequent effect on the stress distribution of the stress field.
As mentioned above, the process of the migration of the industrial wastewater is a complex multiphysics coupled problem, such as thermal-hydraulic (TH), thermalchemical (TC), thermal-mechanical (TM), and THM coupling [8][9][10][11][12][13][14]. Based on the energy conservation model, pore fluid mass conservation equation, and the elastic-plastic constitutive model, a more comprehensive multicomponent THM coupling model for unsaturated porous media was established by Seetharam et al. [15,16]. However, the contribution of the rock formation with strong inhomogeneity to the fluid motion was not considered. et al. established a THMC coupling model for clay barriers in high radioactive waste storage [17]. However, the pore fluid flow and solute diffusion characteristics driven by the combined temperature gradient and concentration gradient were not included. For the study of the mathematical model of fluid seepage, many studies have been carried out. Three typical basic models of rock seepage are proposed, i.e., the equivalent continuous medium seepage model, the fracture network discontinuous medium seepage model, and the dual medium seepage model [18][19][20]. Rossen pointed out that the porefissure and interfracture seepage problems can be solved by the way of Darcy's law and an anisotropic permeability tensor and proposed a method to characterize the anisotropic seepage in fractured rock masses using continuous medium theory [21]. Such equivalence model is not applicable when the fractured rock mass cannot be treated as a continuous medium, and the discontinuous medium method should be employed [22]. Hsieh et al. introduced hydrodynamics into fracture-pore media by using the N-S equation and the continuity equation as governing equations and provided a new approach to solve this difficult problem [23].
Most of the current studies do not simultaneously consider the inhomogeneity of rock mass and the coupling effect of multiphysics during the migration of the contaminants. Therefore, the THM coupling model with the consideration of the multimedia, i.e., matrix-fracture-fault, is proposed. The influence of the fault and natural fractures on the pressure distribution and contaminant concentration during migration is firstly analyzed. Secondly, the influence of fault length, width, dip angle, permeability, and inflow temperature of wastewater on contaminant migration is studied.
The contaminants in industrial wastewater of this paper are not specified as this is a theoretical and parametric study. The contaminant can be referred to as the common heavy metal ion such as Pb (plumbum) or Hg (mercury).

The Coupled THM Model
The seepage of contaminants in rock formations is the result of a combination of physical, chemical, biological, and other factors. Therefore, the study of contaminant migration in rock formations requires the establishment of governing equations under the coupling of multifield such as the chemical field, thermal field, stress field, and hydraulic field. Meanwhile, the migration of contaminants involves convection, adsorption, desorption, hydrodynamic dispersion, chemical action, and microbial decomposition, which cause great difficulties for the research.
In this paper, the study focuses on the seepage mechanism of contaminants in rock formations under the coupled THM effect. The flow of industrial wastewater is considered to be unidirectional, and the effects of adsorption, desorption, chemical interaction, and microbial decomposition during the migration of contaminants are not considered.
In addition, to carry out the numerical modeling, the following simplifications and assumptions of the extremely complex process of THM coupling still required are put forward to establish the theoretical model. (1) The rock is a saturated, isotropic, homogeneous, and small deformation porous elastic medium. (2) The phase change is not considered, i.e., the transformation between gas and liquid phases due to temperature is not considered. (3) Only the heat conduction mode is considered in the rock mass, and thermal radiation is not considered. (4) The thermal strain in the rock mass due to temperature is isotropic.

Equivalent Continuous Medium Model Theory
2.1.1. Hydraulic Field. The density of wastewater is defined as a function of temperature and pressure [24] and can be expressed as where η w is the compression coefficient of the fluid, β w is the thermal expansion coefficient of the fluid, ρ 0 is the fluid reference density, p and p 0 are the pore water pressure and initial pore water pressure, respectively, and T and T 0 are the fluid temperature and initial temperature, respectively. The average value of the general ground temperature gradient is 3°C/100 m; thus, the temperature of the flow of wastewater in the rock formation is below 40°C.
The hydrodynamic viscosity decreases with increasing temperature, and the relationship between the hydrodynamic viscosity and the temperature conforms to the power function. The expression of the dynamic viscosity of 2 Geofluids wastewater with temperature is as follows [25,26]: The flow of industrial wastewater in a rock formation can be assumed as the flow of a fluid in a porous medium. According to the law of conservation of mass and Darcy's law, the hydraulic field equation can be expressed as where c p = 1/ρ∂ρ/∂p,c T = 1/ρ∂ρ/∂T, ϕ is the porosity, ρ is the density of wastewater, k is the permeability of the matrix, g is the acceleration of gravity, ∇H = ð0, 0, 1Þ T , H represents the hydraulic head, and Q m is the source term.
2.1.2. Thermal Field. The fluid flow in the rock mass is assumed in thermal equilibrium. Without the consideration of the heat exchange between the rock mass and the wastewater, the thermal field equation of the fluid flow process in saturated porous media can be obtained as [27] ρC ð Þ eff ∂T ∂t where ðρCÞ eff = ϕρc w + ð1 − ϕÞρ s c s is the effective heat capacity, λ eff = ϕλ w + ð1 − ϕÞλ s is the effective thermal conductivity, α T is the coefficient of thermal expansion of the rock matrix, ε v is the volume strain of the rock matrix, c w is the specific heat capacity of wastewater, λ w is the thermal conductivity of wastewater, ρ s is the density of the rock mass, c s is the specific heat capacity of the rock matrix, and λ s is the thermal conductivity of the rock matrix.

Stress Field.
The constitutive relationship of the rock mass is idealized as isotropic, homogeneous, and linear medium with small deformation. The deformation of the rock mass will be caused by 3 factors: the stress of the rock mass, the pressure of fluid flow in the rock mass, and the temperature change of the rock mass. Therefore, the stressstrain relationship for the rock mass with consideration of pore pressure and temperature is obtained based on the linear principle: where σ ij is the stress tensor, σ ij ′ is the effective stress tensor, G is the shear modulus, G = E/2ð1 + νÞ, ε ij is the strain tensor, λ is the Lame constant,λ = 2Gν/ð1 − 2νÞ, and α is the Biot coefficient: where K is the bulk deformation modulus of the rock matrix, K = E/3ð1 − 2νÞ, and K s is the bulk deformation modulus of the rock grains. Based on the Cauchy strain theory and the static equilibrium, the geometry relationship and the equilibrium equation can be obtained: Combined with Equation (5), we can get the equation of the stress field of rock mass under stress, pore pressure, and temperature: where F i is the volume force of the rock mass and u i is the displacement component of the rock mass.

Porosity Model.
Porosity is the key factor of pore pressure; thus, it is necessary to choose a porosity model of the rock mass. According to the derivation of the definition of porosity, the dynamic porosity equation of the rock mass is defined as follows [28]: 2.1.5. Permeability Model. The Kozeny-Carman (KC) equation [29] is a semiempirical formulation widely used in the study of permeability evolution of the porous media. Based on the porosity model, the permeability model can be expressed as 2.2. THM Coupling Relationship. The coupled THM model defined based on the field equations provided above is illustrated in Figure 1. (1) The hydraulic field affects the mechanical field by changing the fluid pressure.
(2) The mechanical field affects the porosity and permeability through the change of volume strain and thus affects the hydraulic field.  With the continuous improvement of the fracture network, the MFF multimedia model is trending toward perfection in terms of modeling [30].

Hydraulic Field.
The equation controlling the seepage in the fracture can be expressed as follows [28]: where d f is the fracture width, S f is the water storage coefficient of the fissure, ∇ T is the tangential gradient operator along the fissure boundary direction, which is the tangential differentiation on the boundary as a function of the x -direction and y-direction, and k f is the fissure permeability.

Thermal Field.
In contrast to the equivalent continuous medium case, when discrete fractures are present in the model, the discrete fractures become the main flow channels for the fluid in the matrix. Therefore, the effect of convective action in the matrix on the thermal field of the model can be neglected. The heat conduction effect is considered to be present in the matrix only, and therefore, the governing equation for the thermal field in the fracture can be expressed as follows: Since the study in this paper is a 2-dimensional diffusion model, the discrete fractures are treated as a geometry model of line. Therefore, the deformation of the fractures is neglected. The relationship between permeability and stress of the fractures can be established based on the borehole test of water pressure measurement [31]: where k f 0 is the permeability of the fracture when σ n = 0, σ n is the normal corresponding force of the fracture surface, α f is the influence coefficient, and it determines the fracture state.

Contaminant Migration
Model. The migration mechanism of contaminants in porous media consists of 3 parts: convection, i.e., the contaminant is driven by the water flow to migrate downstream; diffusion, where the contaminant diffuses from high concentration area to low concentration area under the effect of a concentration gradient; and dispersion, which is caused by the presence of a porous media, and the porous media cause the migration velocity of contaminants to be different in magnitude and direction from the average flow rate. Although the mechanisms of diffusion and mechanical dispersion are different, they are still difficult to distinguish. And they are generally named as hydrodynamic dispersion. The mechanism of contaminant migration in porous media can be expressed as follows: where c is the solute concentration, D is the hydrodynamic dispersion coefficient, and u is the convection velocity.

Migration in the Multimedia under the Coupled THM Effects
To investigate the migration of the contaminants under the coupled THM effects as well as the influence of the fault and fractures, 2 models are employed, i.e.,

Geofluids
The MF model considers the rock matrix and the fault, and the MFF model considers the rock matrix, the fault, and the fractures.
The upper boundary of the model is the inflow boundary of the wastewater. The initial temperature of the rock formation is 318 K. The temperature of the inflowing contaminant at the upper boundary is 293 K. The concentration of contaminants is 1 mol/m 3 . The head pressure of the upper boundary is 9000 Pa. The parameter used in the numerical study is listed in Table 1.
To study the concentration distribution of the contami-  As shown in Figure 3, the pressure around the fault varies obviously as the contour lines concentrate. The pressure distribution in the model behaves nonhomogeneously due to the fault. The difference in permeability of the fault and the rock matrix may be the reason for this phenomenon as it may cause the nonhomogeneity in the seepage field. Figure 4 is the distribution of the contaminant concentration after 20 years. The contaminants are mainly concentrated in the upper half of the model and along the fault. It may imply that the contaminants migrate very slowly, and the diffusion is more significant in the fault than in the matrix.
As mentioned above, the fault may influence the flow of the wastewater and the migration of the contaminants greatly. The flow vector is employed to illustrate. As shown in Figure 5, both the direction and density of the vectors change obviously around the fault; thus, it can be another aspect to reflect the role of the fault for the distribution of the nonhomogeneity.
Although the fault has a dramatic effect on the flow field, the concentration of the contaminants is monitored along the fault direction, i.e., the line of AB. It is important to note that the length of AB in this paper is calculated based on the start point of A. As shown in Figure 6, the contaminants reach saturation as the concentration is close to 1 mol/m 3 within the length of 3 m toward the fault direction. And   Meanwhile, it can be found that this is a peak in each line of the plot. And the peak occurs at the point of the projection of the endpoint of the fault. The reason may be as follows. The fault has a significant effect on the migration of the contaminants as its permeability is much higher than that of the rock matrix. The projection of the endpoint of the fault is the point which is the nearest to the fault; therefore, the peak value of the concentration occurs around such a point.

MFF Model.
In this section, the MFF model (Figure 8) with the consideration of the rock matrix, the fault, and the fractures is employed to study the effect of the fractures. 50 fractures are generated randomly by using MATLAB coding [32]. The length of the fractures varies from 0.8 m to 4 m. The direction of the fractures varies from 0°to 360°. The aperture of each fracture is 0.0005 m. Compared with Figure 3, there are more drastic pressure changes in the model (Figure 9). The nonhomogeneity of the pressure distribution of the MFF model is more obvious than that of the MF model. The pressure contours change around the natural fractures. However, the main pattern of the pressure distribution of the MFF model is similar to that of the MF model. It can be also found that the distribution of the concentration of the MFF model ( Figure 10) is similar to that of the MF model. The phenomena mentioned above may indicate that the fault has a dominant role in the concentration distribution and contaminant migration.
Based on Figure 11, it can be found that there are also 3 typical zones of the contaminant concentration, the saturation zone, the rapid diffusion zone, and the concentration decrease zone. And the geometry scale of the 3 zones is also similar to that of the MF model. The concentration of the point B reaches 1:15 × 10 −3 mol/m 3 , 1:47 × 10 −3 mol/m 3 , 1:79 × 10 −3 mol/m 3 , 0.016 mol/m 3 , 0.086 mol/m 3 , and    Figure 12 is the plot of the concentration distribution along the bottom of the MFF model (the line of CD); there is also a peak for each line of the plot. Each peak value of MFF is higher than the corresponding one of the MF model. And the position where the peak occurs is almost the same as that of the MFF model. Therefore, the concentration dis-tribution of the 2 models is similar, and the magnitude of the MFF model is higher than that of the MF model.

Effect of the Fault Parameters on Contaminant Migration
As discussed in Section 3, the fault has a dominant role in contaminant migration. The parametric study of the fault  7 Geofluids is carried out in this section. The length, width, dip, permeability, and width of the fault are selected as the parameters. In addition, the effect of the initial temperature of industrial wastewater is also studied.

Length of the Fault.
Based on the results of Section 3.1, the length of the fault may dominate the distribution of contaminants; it is thus a crucial factor affecting the migration of contaminants. In this section, 3 models with different fault lengths (Table 2) are employed to explore the effect of fault length on the migration of contaminants. Figure 13 shows

Fault Dip Angle.
In the process of contaminant migration, the dip angle of the fault may determine the length of the fault; therefore, the migration distance of the contaminants as well as the distribution of the contaminants could be affected by such dip angle. Three models (Table 4) with different dip angles of the fault are considered for the accumulation of contaminants. Figure 15 is        4.5. Inflow Temperature of Industrial Wastewater. When wastewater flows into a rock formation, the inflow temperature will impact the thermal field of the rock formation, and the intereffect of the THM coupling will thus influence the migration of the wastewater. In this section, the inflow temperature of wastewater (Table 6) is selected as a parameter to study the contaminant accumulation.
The accumulation of contaminants along the line of CD under different inflow temperatures is plotted in Figure 17. The accumulation of contaminants at the CD increases gradually with time. After 30 years, the accumulation of contaminants along the line of CD reaches 1.58 mol (Model T1), 1.69 mol (Model T2), and 1.76 mol (Model T3), respectively. The temperature of the inflow water has a positive relationship with contaminant accumulation, while it does not affect as much as the fault parameters. 4.6. Discussion. From Figures 13-17, it can be inferred that the fault length has the most remarkable effect on the contaminant accumulation, while the temperature of the wastewater has the minimum effect on the contaminant accumulation.
Meanwhile, the accumulation of contaminants is divided into 2 stages, slow growth stage and rapid growth stage. And for all the models, the slow growth stage is before 20 years, and the rapid growth stage begins after 20 years. The reason for the 2 stages may be that the main channel of contaminant migration in the slow growth stage is the fault; thus, the concentration accumulation increases slowly. During the rapid growth stage, the long-time scale seepage results in the entire penetration of the contaminants through the    10 Geofluids rock matrix as well as the fault; therefore, the concentration accumulation increases rapidly.

Conclusions
Based on the governing equations of the THM coupled theory, multimedia of MF and MFF is presented. The migration of the contaminants is analyzed. The parametric study of the fault geometry and the temperature of the wastewater is carried out. The conclusion can be summarized as follows: (1) The pressure distribution is strongly nonhomogeneous around the fault. The fracture quantitively affects the concentration distribution, while the fault dominates the concentration distribution and contaminant migration (2) The migration of the contaminants can be geometrically divided into 3 zones along the direction of the fault, the saturation zone, the rapid diffusion zone, and the concentration decrease zone