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.
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 where is the compression coefficient of the fluid, is the thermal expansion coefficient of the fluid, is the fluid reference density, are the pore water pressure and initial pore water pressure, respectively, and 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]:
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 ,, is the porosity, is the density of wastewater, is the permeability of the matrix, is the acceleration of gravity, , represents the hydraulic head, and 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] where is the effective heat capacity, is the effective thermal conductivity, is the coefficient of thermal expansion of the rock matrix, is the volume strain of the rock matrix, is the specific heat capacity of wastewater, is the thermal conductivity of wastewater, is the density of the rock mass, is the specific heat capacity of the rock matrix, and 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: where is the stress tensor, is the effective stress tensor, is the shear modulus, , is the strain tensor, is the Lame constant,, and is the Biot coefficient: where is the bulk deformation modulus of the rock matrix, , and 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 is the volume force of the rock mass and 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]:
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. (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.
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]: where is the fracture width, is the water storage coefficient of the fissure, is the tangential gradient operator along the fissure boundary direction, which is the tangential differentiation on the boundary as a function of the -direction and -direction, and 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:
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
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 is the permeability of the fracture when , is the normal corresponding force of the fracture surface, 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: where is the solute concentration, is the hydrodynamic dispersion coefficient, and 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 model
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/m3. 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 contaminant around the fault and the accumulation of the concentration, 2 monitoring lines ( and ) are selected. The line of is along the direction of the fault. The line of is along the bottom of the geometry models. The coordinates of endpoints 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 Figure 2, a geometric model with the size 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).
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 . It is important to note that the length of in this paper is calculated based on the start point of . 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 ), 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.
It can be also observed that the concentration increases with time. The concentration of the point of reaches mol/m3, mol/m3, 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.
Figure 7 is the plot of the concentration along the line of , 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.
3.2. 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 reaches mol/m3, mol/m3, 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 12 is the plot of the concentration distribution along the bottom of the MFF model (the line of ); 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.
4. 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 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 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 the accumulation of contaminants along the line of under the different fault lengths. The accumulation of contaminants at the line of increases gradually with time. After 30 years, the accumulation of contaminants at the reaches 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.
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 (Table 3) are established to study the influence of fault width on contaminant accumulation.
The accumulation of the contaminants along the line of under the different fault widths is plotted in Figure 14. After 30 years, the accumulation of contaminants at the line of 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.
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 (Table 4) with different dip angles of the fault are considered for the accumulation of contaminants.
Figure 15 is the plot of the accumulation of contaminants along the line of under the different fault dip angles. After 30 years, the accumulation of contaminants along the line of 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.
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 permeability of the fault is set as 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.
Figure 16 shows the accumulation of contaminants along the line of with time under the different permeability. After 30 years, the accumulation of contaminants at 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.
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 under different inflow temperatures is plotted in Figure 17. The accumulation of contaminants at the increases gradually with time. After 30 years, the accumulation of contaminants along the line of 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 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
Symbols
| : | The fluid compression coefficient |
| : | The thermal expansion coefficient of the fluid |
| : | The fluid reference density (kg/m3) |
| : | The density of wastewater (kg/m3) |
| : | The density of the rock mass (kg/m3) |
| : | The pore pressure (MPa) |
| : | The initial pore pressure (MPa) |
| : | The fluid temperature (°C) |
| : | The initial temperature (°C) |
| : | The porosity |
| : | The initial porosity |
| : | The permeability of the matrix (m2) |
| : | The initial permeability of matrix (m2) |
| : | The initial permeability of fault (m2) |
| : | The initial permeability of fracture (m2) |
| : | The acceleration of gravity (m/s2) |
| : | The source term |
| : | The coefficient of thermal expansion of the rock matrix |
| : | The thermal conductivity of the rock matrix (W/(m·°C)) |
| : | The thermal conductivity of wastewater (W/(m·°C)) |
| : | The specific heat capacity of wastewater (J/(kg·°C)) |
| : | The specific heat capacity of the rock matrix (J/(kg·°C)) |
| : | The bulk deformation modulus of the rock matrix (Pa) |
| : | The bulk deformation modulus of the rock grains (Pa) |
| : | Young’s modulus of the rock matrix (Pa) |
| : | Young’s modulus of the rock grains (Pa) |
| : | Poisson’s ratio of the rock matrix |
| : | The volume strain of the rock matrix |
| : | The width of fault (m) |
| : | The aperture of fracture (m). |
Data Availability
The data used to support the findings of this study are available from the first author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This project is supported by the National Natural Science Foundation of China (Nos. 52078477 and 51827901), Key Laboratory of Deep Earth Science and Engineering (Sichuan University), Ministry of Education (DESE202106 and DESE202004), and Guangdong Provincial Key Laboratory of Deep Earth Sciences and Geothermal Energy Exploitation and Utilization (Shenzhen University) (2020-3). The authors acknowledge the supports of the abovementioned funding agencies.