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

Authors: Yanan Gao; Peng Guo; Zetian Zhang; Minghui Li; Feng Gao
Journal: Geofluids
 (2021)

---

## Abstract

Industrial wastewater may have a long-time effect 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 influence of the fault and the fractures on the pressure distribution and contaminant migration is analyzed. The influence 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 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. (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 significant effect on contaminant accumulation. The temperature of the wastewater has the minimum effect 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.


---

## Body

## 1. 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 poisonous and harmful and has adverse effects on the environment and human life [1–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–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), thermal-chemical (TC), thermal-mechanical (TM), and THM coupling [8–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–20]. Rossen pointed out that the pore-fissure 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).


## 2. 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.


### 2.1. 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
(1)ρ=ρ01+ηwp−p0−βwT−T0=ρ01+ηwΔp−βwΔT,where ηw is the compression coefficient of the fluid, βw is the thermal expansion coefficient of the fluid, ρ0 is the fluid reference density, pandp0 are the pore water pressure and initial pore water pressure, respectively, and TandT0 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 wastewater with temperature is as follows [25, 26]:
(2)μ=1.787×10−3e−0.333+1.962×10−4×T×T.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(3)ϕρcp∂p∂t+∇⋅−ρkμ∇p+ρg∇H=Qm−ϕρcT∂T∂t−ρ∂ϕ∂t,where cp=1/ρ∂ρ/∂p,cT=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,1T, H represents the hydraulic head, and Qm 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]
(4)ρCeff∂T∂t+1−ϕTKαT∂εv∂t=λeff∇2T+ρcwkμ∇p∇T,where ρCeff=ϕρcw+1−ϕρscs 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, cw is the specific heat capacity of wastewater, λw is the thermal conductivity of wastewater, ρs is the density of the rock mass, cs is the specific heat capacity of the rock matrix, and λs is the thermal conductivity of the rock matrix.


#### 2.1.3. 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 stress-strain relationship for the rock mass with consideration of pore pressure and temperature is obtained based on the linear principle:(5)σij=σij′+αpδij=2Gεij+λδijεv−KαTδijT−T0+αpδij,where σij is the stress tensor, σij′ is the effective stress tensor, G is the shear modulus, G=E/21+ν, εij is the strain tensor, λ is the Lame constant,λ=2Gν/1−2ν, and α is the Biot coefficient:
(6)α=1−KKs,where K is the bulk deformation modulus of the rock matrix, K=E/31−2ν, and Ks 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:(7)εij=12ui,j+uj,i,(8)εv=εkk=ε11+ε22+ε33,(9)σij,j+fi=0.Combined with Equation (5), we can get the equation of the stress field of rock mass under stress, pore pressure, and temperature:
(10)Gui,ij+uj,ij+λuj,ji−αTT,i+αp,i+Fi=0,where Fi is the volume force of the rock mass and ui is the displacement component of the rock mass.


#### 2.1.4. 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]:
(11)ϕ=α−α−ϕ0exp−1Ksp−p0+αTT−T0−εv−εv0.


#### 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
(12)k=k01ϕ0−1ϕ0−11−1Ksp−p0+αTT−T0−εv−εv03.


### 2.2. THM Coupling Relationship

The coupled THM model defined based on the field equations provided above is illustrated in Figure1. (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. (3) The hydraulic field affects the thermal convection term in the thermal field through the change of fluid velocity vector. (4) The thermal field affects the fluid density by changing the temperature and then affects the hydraulic field. (5) The mechanical field changes the thermal deformation energy to affect the thermal field. (6) The change of temperature in the thermal field produces thermal stress and then affects the mechanical field.Figure 1
THM coupling diagram.


### 2.3. Matrix-Fracture-Fault (MFF) Multimedia Model

The matrix-fracture-fault (MFF) multimedia model is different from the equivalent continuous medium model. It considers the porous medium separately from the fractured medium with discretization of the fractured medium and the influence of each fracture on the hydraulic and thermal fields. With the continuous improvement of the fracture network, the MFF multimedia model is trending toward perfection in terms of modeling [30].


#### 2.3.1. Hydraulic Field

The equation controlling the seepage in the fracture can be expressed as follows [28]:
(13)dfSf∂p∂t+∇T⋅−dfkfμ∇Tp=−df∂εv∂t+Qf,where df is the fracture width, Sf 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 kf is the fissure permeability.


#### 2.3.2. 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:(14)dfρCeff∂T∂t=∇T⋅dfλeff∇TT+dfρcwkμ∇p∇TT+Qw.


#### 2.3.3. Stress Field

Same as the stress governing equation for the equivalent continuous medium, the stress field governing equation for the fracture part can be expressed as(15)Gui,ij+uj,ij+λuj,ji−αTT,i+αp,i+Fi=0.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]:
(16)kf=kf0exp−αfσn,where kf0 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.


### 2.4. 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:(17)∂c∂t=∇⋅D∇c−∇⋅uc+I,where c is the solute concentration, D is the hydrodynamic dispersion coefficient, and u is the convection velocity.


## 2.1. 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
(1)ρ=ρ01+ηwp−p0−βwT−T0=ρ01+ηwΔp−βwΔT,where ηw is the compression coefficient of the fluid, βw is the thermal expansion coefficient of the fluid, ρ0 is the fluid reference density, pandp0 are the pore water pressure and initial pore water pressure, respectively, and TandT0 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 wastewater with temperature is as follows [25, 26]:
(2)μ=1.787×10−3e−0.333+1.962×10−4×T×T.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(3)ϕρcp∂p∂t+∇⋅−ρkμ∇p+ρg∇H=Qm−ϕρcT∂T∂t−ρ∂ϕ∂t,where cp=1/ρ∂ρ/∂p,cT=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,1T, H represents the hydraulic head, and Qm 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]
(4)ρCeff∂T∂t+1−ϕTKαT∂εv∂t=λeff∇2T+ρcwkμ∇p∇T,where ρCeff=ϕρcw+1−ϕρscs 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, cw is the specific heat capacity of wastewater, λw is the thermal conductivity of wastewater, ρs is the density of the rock mass, cs is the specific heat capacity of the rock matrix, and λs is the thermal conductivity of the rock matrix.


### 2.1.3. 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 stress-strain relationship for the rock mass with consideration of pore pressure and temperature is obtained based on the linear principle:(5)σij=σij′+αpδij=2Gεij+λδijεv−KαTδijT−T0+αpδij,where σij is the stress tensor, σij′ is the effective stress tensor, G is the shear modulus, G=E/21+ν, εij is the strain tensor, λ is the Lame constant,λ=2Gν/1−2ν, and α is the Biot coefficient:
(6)α=1−KKs,where K is the bulk deformation modulus of the rock matrix, K=E/31−2ν, and Ks 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:(7)εij=12ui,j+uj,i,(8)εv=εkk=ε11+ε22+ε33,(9)σij,j+fi=0.Combined with Equation (5), we can get the equation of the stress field of rock mass under stress, pore pressure, and temperature:
(10)Gui,ij+uj,ij+λuj,ji−αTT,i+αp,i+Fi=0,where Fi is the volume force of the rock mass and ui is the displacement component of the rock mass.


### 2.1.4. 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]:
(11)ϕ=α−α−ϕ0exp−1Ksp−p0+αTT−T0−εv−εv0.


### 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
(12)k=k01ϕ0−1ϕ0−11−1Ksp−p0+αTT−T0−εv−εv03.


## 2.1.1. Hydraulic Field

The density of wastewater is defined as a function of temperature and pressure [24] and can be expressed as
(1)ρ=ρ01+ηwp−p0−βwT−T0=ρ01+ηwΔp−βwΔT,where ηw is the compression coefficient of the fluid, βw is the thermal expansion coefficient of the fluid, ρ0 is the fluid reference density, pandp0 are the pore water pressure and initial pore water pressure, respectively, and TandT0 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 wastewater with temperature is as follows [25, 26]:
(2)μ=1.787×10−3e−0.333+1.962×10−4×T×T.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(3)ϕρcp∂p∂t+∇⋅−ρkμ∇p+ρg∇H=Qm−ϕρcT∂T∂t−ρ∂ϕ∂t,where cp=1/ρ∂ρ/∂p,cT=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,1T, H represents the hydraulic head, and Qm 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]
(4)ρCeff∂T∂t+1−ϕTKαT∂εv∂t=λeff∇2T+ρcwkμ∇p∇T,where ρCeff=ϕρcw+1−ϕρscs 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, cw is the specific heat capacity of wastewater, λw is the thermal conductivity of wastewater, ρs is the density of the rock mass, cs is the specific heat capacity of the rock matrix, and λs is the thermal conductivity of the rock matrix.


## 2.1.3. 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 stress-strain relationship for the rock mass with consideration of pore pressure and temperature is obtained based on the linear principle:(5)σij=σij′+αpδij=2Gεij+λδijεv−KαTδijT−T0+αpδij,where σij is the stress tensor, σij′ is the effective stress tensor, G is the shear modulus, G=E/21+ν, εij is the strain tensor, λ is the Lame constant,λ=2Gν/1−2ν, and α is the Biot coefficient:
(6)α=1−KKs,where K is the bulk deformation modulus of the rock matrix, K=E/31−2ν, and Ks 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:(7)εij=12ui,j+uj,i,(8)εv=εkk=ε11+ε22+ε33,(9)σij,j+fi=0.Combined with Equation (5), we can get the equation of the stress field of rock mass under stress, pore pressure, and temperature:
(10)Gui,ij+uj,ij+λuj,ji−αTT,i+αp,i+Fi=0,where Fi is the volume force of the rock mass and ui is the displacement component of the rock mass.


## 2.1.4. 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]:
(11)ϕ=α−α−ϕ0exp−1Ksp−p0+αTT−T0−εv−εv0.


## 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
(12)k=k01ϕ0−1ϕ0−11−1Ksp−p0+αTT−T0−εv−εv03.


## 2.2. THM Coupling Relationship

The coupled THM model defined based on the field equations provided above is illustrated in Figure1. (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. (3) The hydraulic field affects the thermal convection term in the thermal field through the change of fluid velocity vector. (4) The thermal field affects the fluid density by changing the temperature and then affects the hydraulic field. (5) The mechanical field changes the thermal deformation energy to affect the thermal field. (6) The change of temperature in the thermal field produces thermal stress and then affects the mechanical field.Figure 1
THM coupling diagram.


## 2.3. Matrix-Fracture-Fault (MFF) Multimedia Model

The matrix-fracture-fault (MFF) multimedia model is different from the equivalent continuous medium model. It considers the porous medium separately from the fractured medium with discretization of the fractured medium and the influence of each fracture on the hydraulic and thermal fields. With the continuous improvement of the fracture network, the MFF multimedia model is trending toward perfection in terms of modeling [30].


### 2.3.1. Hydraulic Field

The equation controlling the seepage in the fracture can be expressed as follows [28]:
(13)dfSf∂p∂t+∇T⋅−dfkfμ∇Tp=−df∂εv∂t+Qf,where df is the fracture width, Sf 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 kf is the fissure permeability.


### 2.3.2. 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:(14)dfρCeff∂T∂t=∇T⋅dfλeff∇TT+dfρcwkμ∇p∇TT+Qw.


### 2.3.3. Stress Field

Same as the stress governing equation for the equivalent continuous medium, the stress field governing equation for the fracture part can be expressed as(15)Gui,ij+uj,ij+λuj,ji−αTT,i+αp,i+Fi=0.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]:
(16)kf=kf0exp−αfσn,where kf0 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.


## 2.3.1. Hydraulic Field

The equation controlling the seepage in the fracture can be expressed as follows [28]:
(13)dfSf∂p∂t+∇T⋅−dfkfμ∇Tp=−df∂εv∂t+Qf,where df is the fracture width, Sf 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 kf is the fissure permeability.


## 2.3.2. 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:(14)dfρCeff∂T∂t=∇T⋅dfλeff∇TT+dfρcwkμ∇p∇TT+Qw.


## 2.3.3. Stress Field

Same as the stress governing equation for the equivalent continuous medium, the stress field governing equation for the fracture part can be expressed as(15)Gui,ij+uj,ij+λuj,ji−αTT,i+αp,i+Fi=0.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]:
(16)kf=kf0exp−αfσn,where kf0 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.


## 2.4. 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:(17)∂c∂t=∇⋅D∇c−∇⋅uc+I,where c is the solute concentration, D is the hydrodynamic dispersion coefficient, and u is the convection velocity.


## 3. 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.,(1)
matrix-fault (MF) model(2)
MFF modelThe 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/m3. The head pressure of the upper boundary is 9000 Pa. The parameter used in the numerical study is listed in Table 1.Table 1
The parameters and values.
ParametersValueYoung’s modulus of the rock matrix,E (Pa)2×1010Young’s modulus of the rock grains,Es (Pa)5×1010Poisson’s ratio,ν0.3Coefficient of thermal expansion of rock matrix,αT (1/K)1.2×10−5Density of the rock grains,ρs (kg/m3)2652Specific heat capacity of the rock matrix,cs (J/(kg·°C))2190Thermal conductivity of the rock matrix,λs (W/(m·°C))0.56Density of the wastewater,ρ0 (kg/m3)1100Specific heat capacity of the wastewater,cw (J/(kg·°C))4210Thermal conductivity of the wastewater,λw (W/(m·°C))0.65Initial porosity,φ00.001Initial permeability of the matrix,k0 (m2)2.5×10−15Initial permeability of the fault,k1 (m2)1×10−13Initial permeability of the fractures,k2 (m2)2.5×10−14Width of the fault,d1 (m)0.002Aperture of the fracture,d2 (m)0.0005To study the concentration distribution of the contaminant around the fault and the accumulation of the concentration, 2 monitoring lines (AB and CD) are selected. The line of AB is along the direction of the fault. The line of CD is along the bottom of the geometry models. The coordinates of endpoints A,B,C,andD are (23.66 m, 30 m), (6.34 m, 0 m), (0 m, 0 m), and (30 m, 0 m).


### 3.1. MF Model

As shown in Figure2, a geometric model with the size 30m×30m is established. The incline angle of the fault is 60°, the width of the fault is 0.002 m, and the coordinate of the midpoint of the fault is (15 m, 15 m).Figure 2
Geometry of the MF model.As shown in Figure3, 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 3
The pressure distribution contour of the MF model.Figure4 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.Figure 4
The concentration distribution contour of the MF model.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 Figure5, 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.Figure 5
The flow vector of the wastewater of the MF model.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 ofAB. 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/m3 within the length of 3 m toward the fault direction. And the diffusion of the contaminants becomes remarkable from about 3 m to 15 m along the fault. Thus, a rapid diffusion zone can be found. From 15 m to the bottom of the model (point B), the concentration decreases rapidly. Therefore, the migration of the contaminants can be divided into 3 zones along the direction of the fault, i.e., the saturation zone, the rapid diffusion zone, and the concentration decrease zone.Figure 6
The contaminant concentration along the line ofAB in the MF model.It can be also observed that the concentration increases with time. The concentration of the point ofB reaches 0.23×10−3 mol/m3, 1.31×10−3 mol/m3, 1.3×10−3 mol/m3, 0.012 mol/m3, 0.074 mol/m3, and 0.138 mol/m3 at the time of 5 years, 10 years, 15 years, 20 years, 25 years, and 30 years, respectively.Figure7 is the plot of the concentration along the line of CD, i.e., the bottom of the model. Similarly, the concentration increases with time, and it becomes the maximum after 30 years. 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.Figure 7
The contaminant concentration along the line ofCD in the MF model.


### 3.2. MFF Model

In this section, the MFF model (Figure8) 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.Figure 8
Geometry of the MFF model.Compared with Figure3, 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.Figure 9
The pressure distribution contour of the MFF model.Figure 10
The concentration distribution contour of the MFF model.Based on Figure11, 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/m3, 1.47×10−3 mol/m3, 1.79×10−3 mol/m3, 0.016 mol/m3, 0.086 mol/m3, and 0.157 mol/m3at the time of 5 years, 10 years, 15 years, 20 years, 25 years, and 30 years, respectively.Figure 11
The contaminant concentration along the line ofAB in the MFF model.Figure12 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 distribution of the 2 models is similar, and the magnitude of the MFF model is higher than that of the MF model.Figure 12
The contaminant concentration along the line ofCD in the MFF model.


## 3.1. MF Model

As shown in Figure2, a geometric model with the size 30m×30m is established. The incline angle of the fault is 60°, the width of the fault is 0.002 m, and the coordinate of the midpoint of the fault is (15 m, 15 m).Figure 2
Geometry of the MF model.As shown in Figure3, 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 3
The pressure distribution contour of the MF model.Figure4 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.Figure 4
The concentration distribution contour of the MF model.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 Figure5, 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.Figure 5
The flow vector of the wastewater of the MF model.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 ofAB. 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/m3 within the length of 3 m toward the fault direction. And the diffusion of the contaminants becomes remarkable from about 3 m to 15 m along the fault. Thus, a rapid diffusion zone can be found. From 15 m to the bottom of the model (point B), the concentration decreases rapidly. Therefore, the migration of the contaminants can be divided into 3 zones along the direction of the fault, i.e., the saturation zone, the rapid diffusion zone, and the concentration decrease zone.Figure 6
The contaminant concentration along the line ofAB in the MF model.It can be also observed that the concentration increases with time. The concentration of the point ofB reaches 0.23×10−3 mol/m3, 1.31×10−3 mol/m3, 1.3×10−3 mol/m3, 0.012 mol/m3, 0.074 mol/m3, and 0.138 mol/m3 at the time of 5 years, 10 years, 15 years, 20 years, 25 years, and 30 years, respectively.Figure7 is the plot of the concentration along the line of CD, i.e., the bottom of the model. Similarly, the concentration increases with time, and it becomes the maximum after 30 years. 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.Figure 7
The contaminant concentration along the line ofCD in the MF model.


## 3.2. MFF Model

In this section, the MFF model (Figure8) 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.Figure 8
Geometry of the MFF model.Compared with Figure3, 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.Figure 9
The pressure distribution contour of the MFF model.Figure 10
The concentration distribution contour of the MFF model.Based on Figure11, 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/m3, 1.47×10−3 mol/m3, 1.79×10−3 mol/m3, 0.016 mol/m3, 0.086 mol/m3, and 0.157 mol/m3at the time of 5 years, 10 years, 15 years, 20 years, 25 years, and 30 years, respectively.Figure 11
The contaminant concentration along the line ofAB in the MFF model.Figure12 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 distribution of the 2 models is similar, and the magnitude of the MFF model is higher than that of the MF model.Figure 12
The contaminant concentration along the line ofCD in the MFF model.


## 4. Effect of the Fault Parameters on Contaminant Migration

As discussed in Section3, the fault has a dominant role in contaminant migration. The parametric study of the fault 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.


### 4.1. Length of the Fault

Based on the results of Section3.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.Table 2
The fault length in the different models.
Model no.Model L1Model L2Model L3Fault length (m)102030Figure13 shows the accumulation of contaminants along the line of CD under the different fault lengths. The accumulation of contaminants at the line of CD increases gradually with time. After 30 years, the accumulation of contaminants at the CD reaches 3.68×10−4 mol (Model L1), 0.064 mol (Model L2), and 1.58 mol (Model L3), respectively. The length of the fault has an obvious effect on the accumulation of contaminants. A longer fault will result in more accumulation of contaminants.Figure 13
The contaminant accumulation along the line ofCD under different fault lengths.


### 4.2. Width of the Fault

The width of a fault may not have a remarkable effect on the migration pattern of the contaminants while it is associated with permeability. It may also have an important effect on contaminant accumulation. In this section, 3 models with the different fault widths (Table3) are established to study the influence of fault width on contaminant accumulation.Table 3
The fault width in the different models.
Model no.Model W1Model W2Model W3Fault width (m)0.0010.0020.003The accumulation of the contaminants along the line ofCD under the different fault widths is plotted in Figure 14. After 30 years, the accumulation of contaminants at the line of CD reaches 0.69 mol (Model W1), 1.58 mol (Model W2), and 1.83 mol (Model W3), respectively. A larger fault width may lead to a larger accumulation. This indicates that the width of the fault can increase the local permeability of the model and enhance the diffusion of contaminants.Figure 14
The contaminant accumulation along the line ofCD under the different fault widths.


### 4.3. 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 (Table4) with different dip angles of the fault are considered for the accumulation of contaminants.Table 4
The fault dip angle in different models.
Model no.Model D1Model D2Model D3Fault dip angle (°)304560Figure15 is the plot of the accumulation of contaminants along the line of CD under the different fault dip angles. After 30 years, the accumulation of contaminants along the line of CD reaches 0.67 mol (Model D1), 1.25 mol (Model D2), and 1.58 mol (Model D3), respectively. The 60° of dip angle of the fault may induce the maximum accumulation of contaminants.Figure 15
The contaminant accumulation along the line ofCD under different fault dip angles.


### 4.4. Fault Permeability

The permeability coefficient is an index of the permeability of rock formations and has a critical influence on the migration of contaminants. In this section, the initial permeabilityk of the fault is set as 1×10−13 m2. Three models (Table 5) with different permeability of the fault are employed to investigate the influence of fault permeability on the accumulation of contaminants.Table 5
Permeability in different models.
Model no.Model P1Model P2Model P3Fault permeability (m2)0.5k1.0k1.5kFigure16 shows the accumulation of contaminants along the line of CD with time under the different permeability. After 30 years, the accumulation of contaminants at CD reaches 0.67 mol (Model P1), 1.58 mol (Model P2), and 1.84 mol (Model P3), respectively. A larger permeability results in a larger accumulation of contaminants.Figure 16
The contaminant accumulation along the line ofCD under different permeability.


### 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 (Table6) is selected as a parameter to study the contaminant accumulation.Table 6
The temperature of the inflow wastewater.
Model no.Model T1Model T2Model T3Inflow temperature (°C)203040The accumulation of contaminants along the line ofCD 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.Figure 17
The contaminant accumulation along the line ofCD under different inflow temperatures.


### 4.6. Discussion

From Figures13–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 rock matrix as well as the fault; therefore, the concentration accumulation increases rapidly.


## 4.1. Length of the Fault

Based on the results of Section3.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.Table 2
The fault length in the different models.
Model no.Model L1Model L2Model L3Fault length (m)102030Figure13 shows the accumulation of contaminants along the line of CD under the different fault lengths. The accumulation of contaminants at the line of CD increases gradually with time. After 30 years, the accumulation of contaminants at the CD reaches 3.68×10−4 mol (Model L1), 0.064 mol (Model L2), and 1.58 mol (Model L3), respectively. The length of the fault has an obvious effect on the accumulation of contaminants. A longer fault will result in more accumulation of contaminants.Figure 13
The contaminant accumulation along the line ofCD under different fault lengths.


## 4.2. Width of the Fault

The width of a fault may not have a remarkable effect on the migration pattern of the contaminants while it is associated with permeability. It may also have an important effect on contaminant accumulation. In this section, 3 models with the different fault widths (Table3) are established to study the influence of fault width on contaminant accumulation.Table 3
The fault width in the different models.
Model no.Model W1Model W2Model W3Fault width (m)0.0010.0020.003The accumulation of the contaminants along the line ofCD under the different fault widths is plotted in Figure 14. After 30 years, the accumulation of contaminants at the line of CD reaches 0.69 mol (Model W1), 1.58 mol (Model W2), and 1.83 mol (Model W3), respectively. A larger fault width may lead to a larger accumulation. This indicates that the width of the fault can increase the local permeability of the model and enhance the diffusion of contaminants.Figure 14
The contaminant accumulation along the line ofCD under the different fault widths.


## 4.3. 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 (Table4) with different dip angles of the fault are considered for the accumulation of contaminants.Table 4
The fault dip angle in different models.
Model no.Model D1Model D2Model D3Fault dip angle (°)304560Figure15 is the plot of the accumulation of contaminants along the line of CD under the different fault dip angles. After 30 years, the accumulation of contaminants along the line of CD reaches 0.67 mol (Model D1), 1.25 mol (Model D2), and 1.58 mol (Model D3), respectively. The 60° of dip angle of the fault may induce the maximum accumulation of contaminants.Figure 15
The contaminant accumulation along the line ofCD under different fault dip angles.


## 4.4. Fault Permeability

The permeability coefficient is an index of the permeability of rock formations and has a critical influence on the migration of contaminants. In this section, the initial permeabilityk of the fault is set as 1×10−13 m2. Three models (Table 5) with different permeability of the fault are employed to investigate the influence of fault permeability on the accumulation of contaminants.Table 5
Permeability in different models.
Model no.Model P1Model P2Model P3Fault permeability (m2)0.5k1.0k1.5kFigure16 shows the accumulation of contaminants along the line of CD with time under the different permeability. After 30 years, the accumulation of contaminants at CD reaches 0.67 mol (Model P1), 1.58 mol (Model P2), and 1.84 mol (Model P3), respectively. A larger permeability results in a larger accumulation of contaminants.Figure 16
The contaminant accumulation along the line ofCD under different permeability.


## 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 (Table6) is selected as a parameter to study the contaminant accumulation.Table 6
The temperature of the inflow wastewater.
Model no.Model T1Model T2Model T3Inflow temperature (°C)203040The accumulation of contaminants along the line ofCD 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.Figure 17
The contaminant accumulation along the line ofCD under different inflow temperatures.


## 4.6. Discussion

From Figures13–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 rock matrix as well as the fault; therefore, the concentration accumulation increases rapidly.


## 5. 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(3)
There is a peak of the concentration along the bottom of the model. Such peak a occurs at the point of the projection of the endpoint of the fault because this position is the nearest to the fault(4)
The accumulation of contaminants is positively correlated with the fault length, width, dip angle, permeability, and temperature of wastewater. The fault length has the most significant effect on contaminant accumulation. The temperature of the wastewater has the minimum effect 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



---

*Source: 5473719-2021-10-16.xml*