A New Analysis Model for Potential Contamination of a Shallow Aquifer from a Hydraulically-Fractured Shale

Hydraulic fracturing (HF) is widely used in shale gas development, which may cause some heavy metals release from shale formations. These contaminants could transport from the fractured shale reservoirs to shallow aquifers. Thus, it is necessary to assess the impact of pollution in shallow aquifers. In this paper, a new analysis model, considering geological distributions, discrete natural fractures (NFs) and faults, is developed to analyze the migration mechanism of contaminants. Furthermore, the alkali erosion of rock caused by high-pH drilling of fluids, is considered in this paper. The numerical results suggest that both NFs and alkali erosion could reduce the time required for contaminants migrating to aquifers. When NFs and alkali erosion are both considered, the migration time will be shortened by 51 years. Alkali erosion makes the impact of NFs, on the contaminant migration, more significant. The migration time decreases with increasing pH values, while the accumulation is on the opposite side. Compared with pH 12.0, the migration time would be increased by 45 years and 29 years for pH 11.0 and 11.5, respectively. However, the migration time for pH 12.5 and 13.0 were found to be decreased by 82 years and 180 years, respectively. Alkali erosion could increase the rock permeability, and the elevated permeability would further enhance the migration velocity of the contaminants, which might play a major role in assessing the potential contamination of shallow aquifers.


Introduction
With the rapid decline in conventional reserves and increase in energy demand, hydraulic fracturing (HF) has emerged as an attractive technology for extracting unconventional resources from low permeability shale reservoirs [1][2][3].However, concerns have arisen about the potential contamination of shallow aquifers [4,5].There are four kinds of contaminants during shale gas extraction: Formation fluid, flow-back fluid, drilling fluid, and HF fluid.The characteristics of these contaminants are described in Table 1.Some investigations [6,7] have shown that a number of heavy metals release from shale formations during gas extractions.Further, experimental evidence [8] shows that flow-back liquid contains some heavy metals.Among these heavy metals, cobalt and zinc are essential nutrients, but can be toxic in larger amounts or certain forms.Other heavy metals, such as cadmium and lead are highly poisonous.The high-pH drilling fluid has also been widely used in the drilling process, which would cause the alkali erosion of rock [9,10].More importantly, the subsurface contains a large number of randomly-oriented natural fractures (NFs), which may enhance the permeability of the rock, locally.Thus, it is urgent to investigate the migration mechanism of heavy metal contaminants in a heterogeneous deep subsurface, considering the alkali erosion of rock.Well, fractures or faults Yang et al. [11] Li et al. [12] Flow-back Fluid Yes 7.0-8.0

Hydraulic Fracturing Fluid and stratigraphic composition
Well, fractures or faults Wang et al. [13] Huang et al. [14] High-pH Drilling Fluid Yes 11.0-12.0 Diesel or mineral oil with a suspended polymer, synthetic polymer or modified clay Drilling fluid loss may occur during drilling, and then filtrate may invades into porous formation Ghavami et al. [15] Kang et al. [16] HF Fluid Yes 6.5-7.0Water, proppant and chemical additives Well, fractures or faults Birdsell et al. [17] Contaminant transport from the fractured shale to aquifers through three potential pathways: Advective transport through bulk media, preferential flow through fractures, and leakage and transport through abandoned wells.Simulations of HF fluid transport from Marcellus shale reservoir to aquifers have been presented by Myers [18], where the contaminants' potential migration pathways are convective transport and fracture transport.He predicted a high risk for the overlying aquifers contamination within 10 years.Gassiat et al. [19] developed a two-dimensional, single-phase, multispecies, density-dependent numerical model for the migration of contaminants.They pointed out that the contaminants from the shale unit could reach the shallow aquifer in less than a 1000 years at 90% of their initial concentration, in the worst-case scenario.Three scenarios for the HF fluid and methane transport from gas reservoir were conducted by Kissinger et al. [20], to determine under which circumstances contaminants would leak into shallower layers.They found that fluid migration is only possible if a combination of several conservative assumptions is met by a scenario.Simulations about the effect of fault on the contaminant migration were conducted by Wei et al. [21], where the distance from the nearest leakage point to the fault, and the inclination of the fault were considered.They found that the migration plume is greater if the distance is closer.Moreover, the fault raises the greatest impact on the migration of contaminants when the dip angle is 70 • .In addition, Reagan et al. [22] presented two general failure scenarios to investigate the migration of gas and water.The two failure scenarios were: Communication via a connecting fracture or fault, and communication via a deteriorated pre-existing, nearby well.The experience of Birdsell et al. [17] had great significance as a reference, which provides a comprehensive perspective on the migration of the HF fluid.However, in the above numerical simulations, the overburden was considered to be homogeneous.It was noted that Pfunt et al. [23] evaluated the risk of HF fluid migration toward shallow aquifers, through a geological subsurface.They found that the injected HF fluid would not reach near-surface aquifers, either during the hydraulic fracturing or in the long term.All these numerical simulations studied the mechanism of contaminant migration from various aspects (see Table 2).Compared with the aforementioned models, this model has the following advantages.Firstly, the overburden is considered to be heterogeneous, which can better reflect the actual geological conditions.Secondly, the existence of NFs is considered in this model, and the interactions of the flow regimes in the matrix, the NFs and the faults is presented for the heterogeneous overburden.Thirdly, the alkali erosion of rock caused by high-pH drilling fluids is considered in this model.

Authors Fluid Driving Force Advantages Disadvantages
Myers [18] HF ∆h imposed at boundaries Taking the lead in studying the migration of HF fluid and providing ideas for later scholars.
Without consideration of the buoyancy caused by difference in density flow.
Gassiat et al. [19] HF Reservoir Overpressure A wide range of parameter studies is presented.
Without consideration of well production.
(1) and (2) do not consider the effects of imbibition and well production; (3) cannot be used for a quantitative risk analysis.
Wei YQ et al. [21] HF Buoyancy The migration of HF caused by density difference is studied in detail.
Without consideration of temperature, alkali erosion and different parameters sensitivity.
Reagan et al. [22] Gas Buoyancy and well The effects of multiphase flow, buoyancy and well production on the migration of HF fluid are considered.

Without consideration of NFs and alkali erosion.
Birdsell et al. [17] HF ∆h imposed at boundaries and well The combination of 5 stages mechanisms is considered and a wide range of parameter studies is presented.
Geological conditions, NFs and alkali erosion are not considered.

Pfunt et al. [23] HF ∆p from injection
The migration of HF caused by density difference is studied in detail.
Without consideration of temperature, alkali erosion and different parameters sensitivity.
In order to evaluate the impact of pollution in shallow aquifers, a new analysis model has been proposed in this paper.Geological distributions, NFs and faults are taken into account to simulate the contaminant migration with a consideration of alkali erosion.Migration mechanism of contaminants in a heterogeneous subsurface is shown in Figure 1.The outline of this paper is as follows.Section 2 discusses the governing equations for contaminant migration in heterogeneous shale gas reservoirs.Section 3 verifies the validity of this model by comparing with an analytical solution, and four simplified cases are given and compared with each other, additionally.Section 4 investigates the effects of NFs and alkali erosion on the migration of contaminants.The findings and conclusions are presented in Section 5.In order to evaluate the impact of pollution in shallow aquifers, a new analysis model has been proposed in this paper.Geological distributions, NFs and faults are taken into account to simulate the contaminant migration with a consideration of alkali erosion.Migration mechanism of contaminants in a heterogeneous subsurface is shown in Figure 1.The outline of this paper is as follows.Section 2 discusses the governing equations for contaminant migration in heterogeneous shale gas reservoirs.Section 3 verifies the validity of this model by comparing with an analytical solution, and four simplified cases are given and compared with each other, additionally.Section 4 investigates the effects of NFs and alkali erosion on the migration of contaminants.The findings and conclusions are presented in Section 5.

Mathematical Model
In this section, the mathematical model for the solute transport is established, which includes groundwater flow and the contaminants' migration.Specially, the interactions of flow regimes in the matrix, the NFs and the faults are presented for a heterogeneous subsurface.In addition, the high-pH drilling fluid would cause the alkali erosion of rock, thus, changing the porosity and permeability of the rock.Therefore, this section mainly introduces the mechanism of contaminant migration in a heterogeneous subsurface, from three aspects: Groundwater flow mechanism, the contaminants' migration mechanism, and the alkali erosion mechanism of shale by a high-pH drilling fluid.Of note, the adsorption, desorption, and microbial decomposition of contaminants have not been considered here.

Groundwater Flow Mechanism
Fluid flow in porous media is very important in a wide range of science and engineering applications, such as shale gas extraction and contaminant migration [24].The mass conservation equation is often used to describe the flow regimes in a porous media.
where φ is the porosity of the porous media, and is a dimensionless number; ρ is the fluid density (kg/m 3 ); v is the seepage velocity (m/s); Q m denotes source/sink term, which refers to the fluid mass that is increased or decreased by external factors, per unit volume per unit time.Darcy's law can be used to describe the flow regime in a shale matrix, and the flow regime in the discrete NFs, follows the cubic law.Compared to the NFs, the traditional Darcy's law is no longer suitable to faults because the groundwater flow in the faults is a nonlinear high-speed flow.Therefore, an equivalent method (equivalent Darcy's law) put forward by Xu [25] is used in this paper to describe the flow regime in faults.The equivalent permeability in faults can be obtained as follows [25]: where K x , K y are the principal permeability of the shale mass with a horizontal fracture (m 2 ), d is the fracture aperture (m); h is the matrix thickness (m); l y = d + 2h; K m is the matrix permeability (m 2 ).When the inclination angle of the fracture is considered, the permeability tensor can be expressed as:

Contaminant Migration Mechanism
The migration mechanism of the contaminants in a porous media include three aspects [26]: Convective, molecular diffusion, and mechanical dispersion.Molecular diffusion and mechanical dispersion are collectively referred to as hydrodynamic dispersion.Convection is the main driving force for the migration of contaminants, and the hydrodynamic dispersion follows Fick's law.The migration mechanism of contaminants, in matrix and fractures, follows the convection-diffusion equation [27].
Energies 2018, 11, 3010 The first term represents the accumulation of the solute concentration, as a function of time.The second and third term represent the amount of change in solute concentration caused by diffusion and convection, respectively.The last term I represents the source/sink term.C is the molar concentration of the solute (mol/m 3 ); t is time (s); D is the hydrodynamic dispersion coefficient (m 2 /s); u is the convection speed (m/s), which is the ratio of the seepage velocity to porosity.
Therefore, Equations ( 1) and ( 5) constitute the basic governing equations for the contaminants' migration in heterogeneous shale gas reservoirs, and they are interdependent.Changes in fluid concentration may cause changes of the fluid density and viscosity, which in turn affects the flow field.

Alkali Erosion Mechanism of Rock by a High-pH Drilling Fluid
The alkali erosion of rock by the high-pH drilling fluid would cause changes in the porosity and the effective bearing surface of the rock [28].SiO 2 is the most abundant substance in the rock that can react with the high-pH drilling fluid.Thus, SiO 2 and OH − are regarded as the reactants in this paper.The reaction equation can be expressed as follows: As the solubility of SiO 2 , in solution, is tiny and basically remains the same, the rate equation for the high-pH drilling fluid can be expressed as: where A is OH − ; B is SiO 2 ; k c denotes reaction rate constant.The relationship between the concentration of OH − (C A ) and reaction rate (V A ) can be expressed as: .
C A is the derivative of C A with respect to time.Assuming that the reaction is carried out at a constant temperature, the expression of the solute concentration over time can be obtained by integrating Equation (8).
where C A0 is the initial concentration of OH − ; t is the duration for which the rock is immersed in the chemical solution.
In addition, the amount of OH − that is consumed by the reactions can be obtained as follows: where V is the volume of the soaking solution.The amount of SiO 2 that is consumed by the reactions can be obtained from Equation (6).
Assuming that the rock components are evenly distributed, the mass of SiO 2 consumed can be expressed as ∆m = N B M B , where M B is the molar mass of SiO 2 .
When the initial mass of the rock is m 0 , the damage variable of the rock, under alkali erosion, can be obtained as follows: Energies 2018, 11, 3010 6 of 22 The evolution of rock porosity φ from its initial value, due to alkali erosion is φ c = φ − φ 0 [29,30].The relationship between the porosity caused by alkali erosion φ c and the damage parameter D c is presented as follows [31]: In addition, the Blake-Kozeny equation can be used to describe the relationship between the permeability and the porosity of the rock [32]: where K is the permeability of the matrix/fracture, K 0 is the initial permeability of the matrix/fracture.Therefore, when the alkali erosion of rock, by high-pH drilling fluids, is considered, the basic governing equations for the contaminant migration can be obtained by Equations ( 1), ( 5), ( 13) and (14).In conclusion, the alkali erosion could increase the rock porosity and permeability, and the elevated rock porosity and permeability would further enhance the velocity of fluid, which is one of the most critical parameters in the calculation of transport equations.

Model Verification
In this section, a new heterogeneous shale gas reservoir model, considering the faults, stratigraphic division, the NFs, and the alkali erosion of rock, is established.The validity of this model was verified through four simplified cases.The four cases, outlined below, were: (1) A single fracture was considered (Case 1); (2) stratigraphic division (Case 2); (3) NFs were considered (Case 3); (4) alkali erosion of rock was considered (Case 4).Governing equations of all cases in this paper were solved by COMSOL 5.2, a commercial software to directly implement partial differential equations.The distribution of the contaminants' concentration in the fracture was analyzed and was compared to the exact solution obtained by Tang [33].The parameters of the four cases are listed in Table 3.
There was a 0.9 m head drop from the lower boundary to the upper boundary, and 1 mol/m 3 contaminant were assumed to flow from the lower boundary.In addition, line 1, whose endpoints were (2,30) and (14,30), was intercepted to research the concentration distribution of the contaminants, at the upper boundary.

Case 1: Consideration of a Single Fracture
In this part, a 30 × 30 m 2 single fracture model is established (see Figure 2a), which is a vertical cross-sectional schematic.The horizontal axis (x-axis) of the graph represents the formation length, while the vertical axis (z-axis) represents the formation depth.The single fracture was replaced by the equivalent region 1 m × 1 m, and the inclination angle of the fracture was 60 • .The equivalent permeability tensor of this fracture could be obtained from Equation (4).

Case 1: Consideration of a Single Fracture
In this part, a 30 × 30 m 2 single fracture model is established (see Figure 2a), which is a vertical cross-sectional schematic.The horizontal axis (x-axis) of the graph represents the formation length, while the vertical axis (z-axis) represents the formation depth.The single fracture was replaced by the equivalent region 1 m × 1 m, and the inclination angle of the fracture was 60°.The equivalent permeability tensor of this fracture could be obtained from Equation ( 4).The location of the contaminants plume, at year 25, is shown in Figure 2b.It clearly shows that the contaminants migrated up along the fault, which might have been caused by a higher permeability in the fault than that in the matrix.After 25 years, a small amount of contaminants reached the upper boundary.This result is consistent with the conclusion of Birdsell [17].
Figure 3 shows the global and local distribution of the contaminants' velocity vector and the streamlines in Model 1, in which the arrows represent the contaminants flow direction and the streamline density reflects the velocity magnitude.It was obvious that the velocity in the fracture was much higher than that in matrix (Figure 3a), where the proportional processing was adopted.In order to show the contaminants flow direction more vividly, normalized processing (see Figure 3b) was employed in this paper.Figure 3c shows the local distribution of the contaminants' velocity vector and the streamlines.As shown in Figure 3, the streamlines in the fault was dense, which also indicated that the velocity in the fault was higher.
Figure 4 shows the concentration curves of contaminants, at the upper boundary.For different times, the concentration evolution of the contaminants had the same tendency, and the contaminants moved to the upper boundary, with a maximum value of 0.1677 mol/m 3 .
Tang [33] proposed an analytical solution for the problem of contaminant transport along a discrete fracture, in a porous rock matrix.In order to compare with the exact analytical solution, the numerical model in Figure 2a was changed to as in Figure 5a, to illustrate the effectiveness of the above results: The inclination angle of the fracture was 90°, the fracture aperture, b = 0.2 cm, the matrix thickness, d = 21.8 cm, and the fracture length, L = 100 cm.The groundwater flow velocity in the fracture was assumed to be 1.0 cm/day, and a contaminant source of 1 mol/m 3 existed at the origin of the fracture.
It can be seen from the Figure 5b that the numerical solution of the solute transport problem is basically consistent with its analytical solution.In general, there is no analytical solution for the The location of the contaminants plume, at year 25, is shown in Figure 2b.It clearly shows that the contaminants migrated up along the fault, which might have been caused by a higher permeability in the fault than that in the matrix.After 25 years, a small amount of contaminants reached the upper boundary.This result is consistent with the conclusion of Birdsell [17].
Figure 3 shows the global and local distribution of the contaminants' velocity vector and the streamlines in Model 1, in which the arrows represent the contaminants flow direction and the streamline density reflects the velocity magnitude.It was obvious that the velocity in the fracture was much higher than that in matrix (Figure 3a), where the proportional processing was adopted.In order to show the contaminants flow direction more vividly, normalized processing (see Figure 3b) was employed in this paper.Figure 3c shows the local distribution of the contaminants' velocity vector and the streamlines.As shown in Figure 3, the streamlines in the fault was dense, which also indicated that the velocity in the fault was higher.
Figure 4 shows the concentration curves of contaminants, at the upper boundary.For different times, the concentration evolution of the contaminants had the same tendency, and the contaminants moved to the upper boundary, with a maximum value of 0.1677 mol/m 3 .
Tang [33] proposed an analytical solution for the problem of contaminant transport along a discrete fracture, in a porous rock matrix.In order to compare with the exact analytical solution, the numerical model in Figure 2a was changed to as in Figure 5a, to illustrate the effectiveness of the above results: The inclination angle of the fracture was 90 • , the fracture aperture, b = 0.2 cm, the matrix thickness, d = 21.8 cm, and the fracture length, L = 100 cm.The groundwater flow velocity in the fracture was assumed to be 1.0 cm/day, and a contaminant source of 1 mol/m 3 existed at the origin of the fracture.
It can be seen from the Figure 5b that the numerical solution of the solute transport problem is basically consistent with its analytical solution.In general, there is no analytical solution for the solute transport, in a fractured rock mass with complex boundaries.Therefore, the mathematical model for the solute transport has a great advantage compared to the analytical method.This model is applicable to the solute transport problems, in fractures, with any inclination angle.It could also be further used to explore the effects of NFs and alkali erosion on contaminant migration.solute transport, in a fractured rock mass with complex boundaries.Therefore, the mathematical model for the solute transport has a great advantage compared to the analytical method.This model is applicable to the solute transport problems, in fractures, with any inclination angle.It could also be further used to explore the effects of NFs and alkali erosion on contaminant migration.solute transport, in a fractured rock mass with complex boundaries.Therefore, the mathematical model for the solute transport has a great advantage compared to the analytical method.This model is applicable to the solute transport problems, in fractures, with any inclination angle.It could also be further used to explore the effects of NFs and alkali erosion on contaminant migration.

Case 2: Consideration of the Geological Distributions
In this part, a simple heterogeneous model is presented, based on Model 1 (see Figure 6a).The formation is divided into two layers, each with a thickness of 15 m.The parameters of each formation could be obtained from Table 2.

Case 2: Consideration of the Geological Distributions
In this part, a simple heterogeneous model is presented, based on Model 1 (see Figure 6a).The formation is divided into two layers, each with a thickness of 15 m.The parameters of each formation could be obtained from Table 2.

Case 2: Consideration of the Geological Distributions
In this part, a simple heterogeneous model is presented, based on Model 1 (see Figure 6a).The formation is divided into two layers, each with a thickness of 15 m.The parameters of each formation could be obtained from Table 2. Figure 7a shows the local distribution of the contaminants' velocity vector and the streamlines in Model 2. Compared with Model 1 (Figure 3c), the streamlines at the fault were relatively sparse, which shows that the flow velocity of the contaminants at the fault, in Model 2, was smaller than that in Model 1. Figure 7b represents the concentration curves of the contaminants at the upper boundary.It was observed that the contaminants moved to the upper boundary, with a maximum value of 0.1518 mol/m 3 , which was smaller than that in Model 1.The above results showed that the Figure 7a shows the local distribution of the contaminants' velocity vector and the streamlines in Model 2. Compared with Model 1 (Figure 3c), the streamlines at the fault were relatively sparse, which shows that the flow velocity of the contaminants at the fault, in Model 2, was smaller than that in Model 1. Figure 7b represents the concentration curves of the contaminants at the upper boundary.It was observed that the contaminants moved to the upper boundary, with a maximum value of 0.1518 mol/m 3 , which was smaller than that in Model 1.The above results showed that the heterogeneity of overburden would affect the migration of the contaminants, which is consistent with that studied by Birdsell [17].
Energies 2018, 11, x FOR PEER REVIEW 10 of 22 heterogeneity of overburden would affect the migration of the contaminants, which is consistent with that studied by Birdsell [17].

Case 3: Consideration of the NFs
In this subsection, 45 NFs were randomly embedded into the matrix, based on Model 2 and the NFs' aperture, d f was 1.0 × 10 −4 m (see Figure 8a).These fractures were simplified into one-dimensional line element, which could be created by MATLAB (2015).

Case 3: Consideration of the NFs
In this subsection, 45 NFs were randomly embedded into the matrix, based on Model 2 and the NFs' aperture, f d was 1.0 × 10 −4 m (see Figure 8a).These fractures were simplified into one-dimensional line element, which could be created by MATLAB (2015).Figure 8b shows the location of the contaminants' plume, at year 25, in Model 3. Compared with Model 2, more contaminants migrated to the upper boundary.The contaminants concentration at the bottom of the fault in Model 3 was less than that in Model 2. These results showed that NFs could increase the local permeability of the matrix, provided more contaminants migrated to the upper boundary.
Figure 9a shows the local distribution of the contaminants' velocity vector in Model 3. As shown in Figure 9a, the contaminants' flowed from the matrix to the NFs and the fault, because there was a higher permeability, which is consistent with that studied by Zhang [34]. Figure 9b represents the concentration curves of the contaminants at the upper boundary.It was observed that the  Figure 9a shows the local distribution of the contaminants' velocity vector in Model 3. As shown in Figure 9a, the contaminants' flowed from the matrix to the NFs and the fault, because there was a higher permeability, which is consistent with that studied by Zhang [34]. Figure 9b represents the concentration curves of the contaminants at the upper boundary.It was observed that the contaminants moved to the upper boundary with a maximum value of 0.2128 mol/m 3 , which was higher than that in Model 1 and Model 2. This shows that the NFs had a significant effect on the amount of contaminants reaching the upper boundary.OH− , which could be obtained from Kang et al. [16].In addition, the increment of rock porosity with time caused by the alkali erosion (φ c ) could be fitted with a third order polynomial, to facilitate the calculation by COMSOL 5.2.
Figure 10a shows the increment of porosity and its fitting curve.As shown in Figure 10a, the porosity of the rock increased by 0.00514, about 150 days later, and the fitting curve could well reflect the change in the rock porosity, caused by erosion.In addition, assuming that the matrix's initial permeability was 5.0 × 10 −16 m 2 , the increment of permeability is shown in Figure 10b

Case 4: Alkali Erosion of Rock by High-pH Drilling Fluids
Some high-pH (pH = 11-12) drilling fluids are widely used in the drilling process to prevent the wellbore instability, which could cause an alkali erosion of the rock.pH = 12 was chosen in this section to study the influence of alkali erosion on the contaminants' migration.When pH = 12.0, the initial concentration of OH − was 0.01 mol/L.The relationship between the average reaction rate   Figure 10a shows the increment of porosity and its fitting curve.As shown in Figure 10a, the porosity of the rock increased by 0.00514, about 150 days later, and the fitting curve could well reflect the change in the rock porosity, caused by erosion.In addition, assuming that the matrix's initial permeability was 5.0 × 10 −16 m 2 , the increment of permeability is shown in Figure 10b Figure 11a shows the location of the contaminants' plume, at year 25, in Model 4. Compared to Model 3, more contaminants migrated to the top of the fault and the upper boundary.In addition, there were more contaminants around the NFs. Figure 11b represents the concentration curves of contaminants at the upper boundary.It was observed that the contaminants moved to the upper boundary with a maximum value of 0.2371 mol/m 3 , which was higher than that in Model 3.This showed that the alkali erosion had a significant effect on the amount of contaminants reaching the upper boundary.
Model 3, more contaminants migrated to the top of the fault and the upper boundary.In addition, there were more contaminants around the NFs. Figure 11b represents the concentration curves of contaminants at the upper boundary.It was observed that the contaminants moved to the upper boundary with a maximum value of 0.2371 mol/m 3 , which was higher than that in Model 3.This showed that the alkali erosion had a significant effect on the amount of contaminants reaching the upper boundary.Figure 12 shows the comparison of the contaminants concentration in the four cases.As shown in Figure 12a, the concentration of contaminants at point (15−5 3 , 30) was 0.09 mol/m 3 , 0.07 mol/m 3 , 0.11 mol/m 3 and 0.14 mol/m 3 for the four models, after 30 years, respectively.Figure 12b represents the concentration curves of the contaminants at the upper boundary for the four models, after 30 years.The concentration of the contaminants, at the same position, was highest in Model 4, followed by Model 3, Model 1, and, finally, Model 2. This showed that the NFs and the alkali erosion both had a great influence on the migration of contaminants.Figure 12 shows the comparison of the contaminants concentration in the four cases.As shown in Figure 12a, the concentration of contaminants at point (15−5 √ 3, 30) was 0.09 mol/m 3 , 0.07 mol/m 3 , 0.11 mol/m 3 and 0.14 mol/m 3 for the four models, after 30 years, respectively.Figure 12b represents the concentration curves of the contaminants at the upper boundary for the four models, after 30 years.The concentration of the contaminants, at the same position, was highest in Model 4, followed by Model 3, Model 1, and, finally, Model 2. This showed that the NFs and the alkali erosion both had a great influence on the migration of contaminants.

Model Application
Model 3, more contaminants migrated to the top of the fault and the upper boundary.In addition, there were more contaminants around the NFs. Figure 11b represents the concentration curves of contaminants at the upper boundary.It was observed that the contaminants moved to the upper boundary with a maximum value of 0.2371 mol/m 3 , which was higher than that in Model 3.This showed that the alkali erosion had a significant effect on the amount of contaminants reaching the upper boundary.Figure 12 shows the comparison of the contaminants concentration in the four cases.As shown in Figure 12a, the concentration of contaminants at point (15−5 3 , 30) was 0.09 mol/m 3 , 0.07 mol/m 3 , 0.11 mol/m 3 and 0.14 mol/m 3 for the four models, after 30 years, respectively.Figure 12b represents the concentration curves of the contaminants at the upper boundary for the four models, after 30 years.The concentration of the contaminants, at the same position, was highest in Model 4, followed by Model 3, Model 1, and, finally, Model 2. This showed that the NFs and the alkali erosion both had a great influence on the migration of contaminants.

Model Application
In this section, a new analysis model, considering the heterogeneity of the subsurface and the alkali erosion of rock, was chosen as an application example to analyze the effect of the NFs and the alkali erosion, on the aquifer's pollution levels.Darcy's law, the cubic law, and the equivalent Darcy's law were applied in the matrix, the discrete NFs and the faults, respectively.In addition, three kinds of excessive heavy metals (Cd, Pb and Co) which come from the Cambrian Niutitang Formation flow-back liquid, were selected as the contaminants.Information about the three contaminants is shown in Table 4.

Geometric Model
As Figure 13a shows, an 8000 × 2800 m 2 deep subsurface model, containing 500 discrete NFs and 7 hydraulic fractures, was established.Here the length of the fractures was about 20 m-300 m.
Fractures were simplified into one-dimensional line element, and randomly distributed natural fracture network was then created using MATLAB (2015).The geometric model was meshed by a random triangle mesh, and the local grid refinement was employed in the fracture region, as shown in Figure 13b.An 84 m head drop was assumed from the lower boundary to the upper boundary.In addition, kinds of contaminants were assumed to flow from the lower boundary, they were 689.09ug/L Cd, 10.28 ug/L Pb, and 200.15 ug/L Co, respectively.
Energies 2018, 11, x FOR PEER REVIEW 13 of 22 In this section, a new analysis model, considering the heterogeneity of the subsurface and the alkali erosion of rock, was chosen as an application example to analyze the effect of the NFs and the alkali erosion, on the aquifer's pollution levels.Darcy's law, the cubic law, and the equivalent Darcy's law were applied in the matrix, the discrete NFs and the faults, respectively.In addition, three kinds of excessive heavy metals (Cd, Pb and Co) which come from the Cambrian Niutitang Formation flow-back liquid, were selected as the contaminants.Information about the three contaminants is shown in Table 4.

Geometric Model
As Figure 13a shows, an 8000 × 2800 m 2 deep subsurface model, containing 500 discrete NFs and 7 hydraulic fractures, was established.Here the length of the fractures was about 20 m-300 m.
Fractures were simplified into one-dimensional line element, and randomly distributed natural fracture network was then created using MATLAB (2015).The geometric model was meshed by a random triangle mesh, and the local grid refinement was employed in the fracture region, as shown in Figure 13b.An 84 m head drop was assumed from the lower boundary to the upper boundary.In addition, three kinds of contaminants were assumed to flow from the lower boundary, they were 689.09ug/L Cd, 10.28 ug/L Pb, and 200.15 ug/L Co, respectively.

Model Parameters and Properties
The permeability (K) and the effective porosity for each of the stratigraphic units are shown in Table 5.These were obtained from the existing literature and adjusted, but were still consistent with the realistic permeability of the overburden.The general parameters of this model are listed in Table 6.

Numerical Results
Two useful metrics were employed in this section to track the contaminant migration: (i) The concentration of the contaminants at the top of the fault, which was useful for identifying at which point of time the contaminants came into contact with the aquifers; (ii) the accumulation of the heavy metal contaminants in the aquifers, which was used to evaluate the water quality of the aquifers.
Figure 14a shows the concentration curves of the three kinds of contaminants, with the time-point at the top of the fault.As Figure 14a shows, the three kinds of contaminants migrated to the aquifers about after 231 years.After a 1000 years, the concentration of Cd, Pb, and Co would increase to 0.005 mol/m 3 , 4.027 × 10 −5 mol/m 3 , and 0.00281 mol/m 3 , respectively.Figure 14b shows the variation curves of the heavy metal contaminants accumulation, with time, in the aquifers.After a 1000 years, the accumulation of Cd, Pb, and Co, in the aquifers would increase to 106.48 mol, 0.8568 mol, and 60.84 mol, respectively.
the aquifers about after 231 years.After a 1000 years, the concentration of Cd, Pb, and Co would increase to 0.005 mol/m 3 , 4.027 × 10 −5 mol/m 3 , and 0.00281 mol/m 3 , respectively.Figure 14b shows the variation curves of the heavy metal contaminants accumulation, with time, in the aquifers.After a 1000 years, the accumulation of Cd, Pb, and Co, in the aquifers would increase to 106.48 mol, 0.8568 mol, and 60.84 mol, respectively.As the migration mechanism of the three contaminants was consistent, Cd was considered, as an example, to study the migration mechanism of the contaminants in the latter studies.
The slice plots in Figure 15 represent the concentration distribution of the contaminants at 100 years, 250 years, 500 years, and 1000 years, respectively.As shown in Figure 15, the contaminants spread along the fault and a small amount of contaminants was found entering the aquifers at 250 years.Subsequently, the contaminants began to migrate to the aquifers and the area of the contaminated aquifers expanded gradually.As the migration mechanism of the three contaminants was consistent, Cd was considered, as an example, to study the migration mechanism of the contaminants in the latter studies.
The slice plots in Figure 15 represent the concentration distribution of the contaminants at 100 years, 250 years, 500 years, and 1000 years, respectively.As shown in Figure 15, the contaminants spread along the fault and a small amount of contaminants was found entering the aquifers at 250 years.Subsequently, the contaminants began to migrate to the aquifers and the area of the contaminated aquifers expanded gradually.
the aquifers about after 231 years.After a 1000 years, the concentration of Cd, Pb, and Co would increase to 0.005 mol/m 3 , 4.027 × 10 −5 mol/m 3 , and 0.00281 mol/m 3 , respectively.Figure 14b shows the variation curves of the heavy metal contaminants accumulation, with time, in the aquifers.After a 1000 years, the accumulation of Cd, Pb, and Co, in the aquifers would increase to 106.48 mol, 0.8568 mol, and 60.84 mol, respectively.As the migration mechanism of the three contaminants was consistent, Cd was considered, as an example, to study the migration mechanism of the contaminants in the latter studies.
The slice plots in Figure 15 represent the concentration distribution of the contaminants at 100 years, 250 years, 500 years, and 1000 years, respectively.As shown in Figure 15, the contaminants spread along the fault and a small amount of contaminants was found entering the aquifers at 250 years.Subsequently, the contaminants began to migrate to the aquifers and the area of the contaminated aquifers expanded gradually.Figure 16 shows the local distribution of the contaminants' velocity vector.As shown in Figure 16a, the contaminants flowed into the fault, along the hydraulic fractures.The contaminants formed a preferential flow in the fractures and the fault because they had a higher permeability, which demonstrated that the heterogeneity of the reservoir had a significant effect on the flow of the contaminants.Line EF, whose endpoints were (4389.5, −2506.7)and (4476.1,−2456.7),was intercepted to discuss the velocity distribution of the contaminants, perpendicular to the fault direction.As shown in Figure 16b, the flow rate of the contaminants in the fault was three orders of magnitude larger than that in the matrix.
The migration mechanisms of the contaminants in the porous media included three aspects: (i) Convective, (ii) molecular diffusion, (iii) mechanical dispersion.Figure 17a shows the variation curves of the Cd flux, with time, caused by these three effects, at the top of the fault.Convection played a dominant role in the migration of the contaminants, and the flux of contaminants, caused by it, was 3-4 orders of magnitude larger than that of the other two.Figure 17b shows the variation in the percentage of Cd, with time, in each region.After a 1000 years, the percentage of the Cd-content in the three regions was 0.879, 0.0703, and 0.0507, respectively, which showed that although most of the contaminants were retained in the subsurface, some contaminants still reached the aquifers.Figure 16 shows the local distribution of the contaminants' velocity vector.As shown in Figure the contaminants flowed into the fault, along the hydraulic fractures.The contaminants formed a preferential flow in the fractures and the fault because they had a higher permeability, which demonstrated that the heterogeneity of the reservoir had a significant effect on the flow of the contaminants.Line EF, whose endpoints were (4389.5, −2506.7)and (4476.1,−2456.7),was intercepted to discuss the velocity distribution of the contaminants, perpendicular to the fault direction.As shown in Figure 16b, the flow rate of the contaminants in the fault was three orders of magnitude larger than that in the matrix.
The migration mechanisms of the contaminants in the porous media included three aspects: (i) Convective, (ii) molecular diffusion, (iii) mechanical dispersion.Figure 17a shows the variation curves of the Cd flux, with time, caused by these three effects, at the top of the fault.Convection played a dominant role in the migration of the contaminants, and the flux of contaminants, caused by it, was 3-4 orders of magnitude larger than that of the other two.Figure 17b shows the variation in the percentage of Cd, with time, in each region.After a 1000 years, the percentage of the Cd-content in the three regions was 0.879, 0.0703, and 0.0507, respectively, which showed that although most of the contaminants were retained in the subsurface, some contaminants still reached the aquifers.Figure 16 shows the local distribution of the contaminants' velocity vector.As shown in Figure 16a, the contaminants flowed into the fault, along the hydraulic fractures.The contaminants formed a preferential flow in the fractures and the fault because they had a higher permeability, which demonstrated that the heterogeneity of the reservoir had a significant effect on the flow of the contaminants.Line EF, whose endpoints were (4389.5, −2506.7)and (4476.1,−2456.7),was intercepted to discuss the velocity distribution of the contaminants, perpendicular to the fault direction.As shown in Figure 16b, the flow rate of the contaminants in the fault was three orders of magnitude larger than that in the matrix.
The migration mechanisms of the contaminants in the porous media included three aspects: (i) Convective, (ii) molecular diffusion, (iii) mechanical dispersion.Figure 17a shows the variation curves of the Cd flux, with time, caused by these three effects, at the top of the fault.Convection played a dominant role in the migration of the contaminants, and the flux of contaminants, caused by it, was 3-4 orders of magnitude larger than that of the other two.Figure 17b shows the variation in the percentage of Cd, with time, in each region.After a 1000 years, the percentage of the Cd-content in the three regions was 0.879, 0.0703, and 0.0507, respectively, which showed that although most of the contaminants were retained in the subsurface, some contaminants still reached the aquifers.The variation curves of the concentration and accumulation of Cd, with time, for the three cases are shown in Figure 18a,b, respectively.As shown in Figure 18a, the time required for the contaminants' migration to the aquifers was 282 years, 268 years, and 231 years, for the three scenarios, respectively.It demonstrated that both the NFs and the alkali erosion had a great influence on the contaminants' migration.The alkali erosion made the impact of the NFs on the contaminants' migration more significant.In addition, the maximum accumulation of Cd in the aquifers was 80.18 mol, 92.68 mol, and 106.48 mol for the three scenarios, respectively.In this section, three cases (see Table 7) are presented, based on the existence of the NFs and the alkali erosion: (1) Neither NFs nor alkali erosion; (2) with NFs, without alkali erosion; (3) both NFs and alkali erosion.The variation curves of the concentration and accumulation of Cd, with time, for the three cases are shown in Figure 18a,b, respectively.As shown in Figure 18a, the time required for the contaminants' migration to the aquifers was 282 years, 268 years, and 231 years, for the three scenarios, respectively.It demonstrated that both the NFs and the alkali erosion had a great influence on the contaminants' migration.The alkali erosion made the impact of the NFs on the contaminants' migration more significant.In addition, the maximum accumulation of Cd in the aquifers was 80.18 mol, 92.68 mol, and 106.48 mol for the three scenarios, respectively.The existence of the NFs could increase the local permeability of the matrix, which had a great influence on the contaminants' migration.In this section, different number of NFs were selected to explore the effect of the NFs' quantity on the migration of the contaminants.They were 300 NFs, 500 NFs and 700 NFs, respectively.
Figure 19a,b show the variation curves of the concentration and accumulation of Cd, with time, for the different number of NFs, respectively.As shown in Figure 19, both the concentration and the accumulation of Cd increased with the number of NFs.The time required for the contaminants migrating to the aquifers was 238 years, 231 years, and 214 years, for 300 NFs, 500 NFs, and 700 NFs, respectively.In addition, the maximum accumulation of Cd in aquifers was 100.43 mol, 106.48 mol, and 113.39 mol for 300 NFs, 500 NFs, and 700 NFs, respectively.The more the number of NFs, the shorter was the time the contaminants took to migrate to the aquifers, and the greater the accumulation in the aquifers.The existence of the NFs could increase the local permeability of the matrix, which had a great influence on the contaminants' migration.In this section, different number of NFs were selected to explore the effect of the NFs' quantity on the migration of the contaminants.They were 300 NFs, 500 NFs and 700 NFs, respectively.
Figure 19a,b show the variation curves of the concentration and accumulation of Cd, with time, for the different number of NFs, respectively.As shown in Figure 19, both the concentration and the accumulation of Cd increased with the number of NFs.The time required for the contaminants migrating to the aquifers was 238 years, 231 years, and 214 years, for 300 NFs, 500 NFs, and 700 NFs, respectively.In addition, the maximum accumulation of Cd in aquifers was 100.43 mol, 106.48 mol, and 113.39 mol for 300 NFs, 500 NFs, and 700 NFs, respectively.The more the number of NFs, the shorter was the time the contaminants took to migrate to the aquifers, and the greater the accumulation in the aquifers.could be obtained from Kang et al. [16].In addition, the general expression of the porosity could be fitted with a simple third order polynomial: The correlation coefficients for the different pH values are shown in Table 8.Note: B1, B2, and B3 denote the correlation coefficients for a third order polynomial fit.B1 is the coefficient of the first term; B2 is the quadratic coefficient; B3 is the cubic coefficient.I here refers to the vertical intercept, which is a point where the graph of a function intersects the y-axis of the coordinate system.As such, the point satisfies x = 0.

The Alkali Erosion on the Water Quality of the Aquifers
In this section, five different pH values (pH = 11.0,11.5, 12.0, 12.5, 13.0) were used to study the effect of alkali erosion on the contaminants' migration.The relationship between the average reaction rate (k c ) and the pH C OH − could be obtained from Kang et al. [16].In addition, the general expression of the porosity could be fitted with a simple third order polynomial: The correlation coefficients for the different pH values are shown in Table 8.Note: B 1 , B 2 , and B 3 denote the correlation coefficients for a third order polynomial fit.B 1 is the coefficient of the first term; B 2 is the quadratic coefficient; B 3 is the cubic coefficient.I here refers to the vertical intercept, which is a point where the graph of a function intersects the y-axis of the coordinate system.As such, the point satisfies x = 0.

Conclusions
Shale gas extraction by hydraulic fracturing has raised concerns about the potential contamination of shallow aquifers.In this paper, a new analysis model has been proposed for the contaminants' migration, in a heterogeneous subsurface.Geological distributions, NFs, and faults were taken into account to simulate the contaminants' migration, with consideration of the alkali erosion.Based on this model, a two-dimensional heterogeneous subsurface was numerically studied for the impact of pollution in shallow aquifers, caused by shale gas extraction.Then, the influences of NFs and alkali erosion on the contaminants migration were investigated through two useful metrics.One of which was the contaminants concentration at the top of the fault; the other was the contaminants accumulation in the aquifer.The present numerical simulations revealed that both the NFs and the alkali erosion could reduce the time required for the contaminants' migration to the

Conclusions
Shale gas extraction by hydraulic fracturing has raised concerns about the potential contamination of shallow aquifers.In this paper, a new analysis model has been proposed for the contaminants' migration, in a heterogeneous subsurface.Geological distributions, NFs, and faults were taken into account to simulate the contaminants' migration, with consideration of the alkali erosion.Based on this model, a two-dimensional heterogeneous subsurface was numerically studied for the impact of pollution in shallow aquifers, caused by shale gas extraction.Then, the influences of NFs and alkali erosion on the contaminants migration were investigated through two useful metrics.One of which was the contaminants concentration at the top of the fault; the other was the contaminants accumulation in the aquifer.The present numerical simulations revealed that both the NFs and the alkali erosion could reduce the time required for the contaminants' migration to the

Conclusions
Shale gas extraction by hydraulic fracturing has raised concerns about the potential contamination of shallow aquifers.In this paper, a new analysis model has been proposed for the contaminants' migration, in a heterogeneous subsurface.Geological distributions, NFs, and faults were taken into account to simulate the contaminants' migration, with consideration of the alkali erosion.Based on this model, a two-dimensional heterogeneous subsurface was numerically studied for the impact of pollution in shallow aquifers, caused by shale gas extraction.Then, the influences of NFs and alkali erosion on the contaminants migration were investigated through two useful metrics.One of which was the contaminants concentration at the top of the fault; the other was the contaminants accumulation in the aquifer.The present numerical simulations revealed that both the NFs and the alkali erosion could reduce the time required for the contaminants' migration to the aquifers.The migration time decreased with an increase in the number of NFs, while the accumulation was on the opposite side.Compared with 500 NFs, the migration time for 700 NFs would be shortened by 17 years, while the migration time for 300 NFs would be increased by 7 years.This information could be used to recognize the influence of natural fractures on contaminants' migration.Alkali erosion of rock, by a high-pH drilling fluid, results in the increase of rock porosity and permeability, which would have a significant effect on the contaminants' migration.The migration time decreased with the increase of the pH values, while the accumulation was on the opposite side.Compared with a pH of 12.0, the migration time would be increased by 45 years and 29 years for pH 11.0 and 11.5, respectively.However, the migration time for pH 12.5 and 13.0 were found to be decreased by 82 years and 180 years, respectively.In particular, when the NFs and the alkali erosion were both considered, the migration time would be shortened by 51 years.Additionally, similar to that shown in previous studies, three heavy metal contaminants migrated up, along the fault.Convection played a dominant role in the process of the contaminants' migration, and the contaminants' flux, caused by it was 3-4 orders of magnitude larger, than that of the hydrodynamic dispersion.In conclusion, the results of the present study showed that alkali erosion of shale could increase rock permeability, which further enhanced the migration velocity of the contaminants.This could provide new insights into the assessment of potential contamination in shallow aquifers.The migration mechanism of contaminants mainly includes convection, hydrodynamic dispersion, adsorption, desorption, microbial decomposition, and the chemical reaction.However, the adsorption, desorption, and microbial decomposition of contaminants were not considered in this paper, which should be studied further.

Figure 1 .
Figure 1.Migration mechanism of contaminants in a heterogeneous subsurface.

Figure 1 .
Figure 1.Migration mechanism of contaminants in a heterogeneous subsurface.

Figure 3 .
Figure 3.The contaminants velocity vector and the streamlines in Model 1.(a) Local distribution of the contaminants' velocity vector, where proportional processing is employed.(b) Global distribution of the contaminants' velocity vector, where normalized processing is employed.(c) Local distribution of the contaminants' velocity vector and the streamlines.

Figure 4 .
Figure 4. Contaminants concentration at the upper boundary, at different times.

Figure 3 .
Figure 3.The contaminants velocity vector and the streamlines in Model 1.(a) Local distribution of the contaminants' velocity vector, where proportional processing is employed.(b) Global distribution of the contaminants' velocity vector, where normalized processing is employed.(c) Local distribution of the contaminants' velocity vector and the streamlines.

Figure 3 .
Figure 3.The contaminants velocity vector and the streamlines in Model 1.(a) Local distribution of the contaminants' velocity vector, where proportional processing is employed.(b) Global distribution of the contaminants' velocity vector, where normalized processing is employed.(c) Local distribution of the contaminants' velocity vector and the streamlines.

Figure 4 .
Figure 4. Contaminants concentration at the upper boundary, at different times.

Figure 4 .Figure 5 .
Figure 4. Contaminants concentration at the upper boundary, at different times.Energies 2018, 11, x FOR PEER REVIEW 9 of 22

Figure 5 .
Figure 5. (a) Fracture-matrix system.(b) Comparison of the analytical solution and the numerical solution.

Figure 5 .
Figure 5. (a) Fracture-matrix system.(b) Comparison of the analytical solution and the numerical solution.

Figure 6 .
Figure 6.(a) Geometric Model 2. (b) Concentration distribution of contaminants in Model 2. The location of the contaminants' plume at year 25, in Model 2, is shown in Figure 6b.Compared to Model 1, more contaminants accumulated at the bottom of the fault in Model 2. This indicated that the contaminants in Model 2 was more inclined to flow to the fault, compared with Model 1, which might have been caused by the lower permeability of the underlying layer in Model 2.Figure7ashows the local distribution of the contaminants' velocity vector and the streamlines in Model 2. Compared with Model 1 (Figure3c), the streamlines at the fault were relatively sparse, which shows that the flow velocity of the contaminants at the fault, in Model 2, was smaller than that in Model 1. Figure7brepresents the concentration curves of the contaminants at the upper boundary.It was observed that the contaminants moved to the upper boundary, with a maximum value of 0.1518 mol/m 3 , which was smaller than that in Model 1.The above results showed that the

Figure 6 .
Figure 6.(a) Geometric Model 2. (b) Concentration distribution of contaminants in Model 2. The location of the contaminants' plume at year 25, in Model 2, is shown in Figure 6b.Compared to Model 1, more contaminants accumulated at the bottom of the fault in Model 2. This indicated that the contaminants in Model 2 was more inclined to flow to the fault, compared with Model 1, which might have been caused by the lower permeability of the underlying layer in Model 2.Figure7ashows the local distribution of the contaminants' velocity vector and the streamlines in Model 2. Compared with Model 1 (Figure3c), the streamlines at the fault were relatively sparse, which shows that the flow velocity of the contaminants at the fault, in Model 2, was smaller than that in Model 1. Figure7brepresents the concentration curves of the contaminants at the upper boundary.It was observed that the contaminants moved to the upper boundary, with a maximum value of 0.1518 mol/m 3 , which was smaller than that in Model 1.The above results showed that the heterogeneity of overburden would affect the migration of the contaminants, which is consistent with that studied by Birdsell[17].

Figure 7 .
Figure 7. (a) Contaminants velocity vector and the streamlines.(b) Contaminants' concentration at the upper boundary, in Model 2.

3. 3 .
Case 3: Consideration of the NFs In this subsection, 45 NFs were randomly embedded into the matrix, based on Model 2 and the NFs' aperture, f d was 1.0 × 10 −4 m (see Figure 8a).These fractures were simplified into one-dimensional line element, which could be created by MATLAB (2015).

Figure 7 .
Figure 7. (a) Contaminants velocity vector and the streamlines.(b) Contaminants' concentration at the upper boundary, in Model 2.

Figure 7 .
Figure 7. (a) Contaminants velocity vector and the streamlines.(b) Contaminants' concentration at the upper boundary, in Model 2.

Figure
Figure8bshows the location of the contaminants' plume, at year 25, in Model 3. Compared with Model 2, more contaminants migrated to the upper boundary.The contaminants concentration at the bottom of the fault in Model 3 was less than that in Model 2. These results showed that NFs could increase the local permeability of the matrix, provided more contaminants migrated to the upper boundary.Figure9ashows the local distribution of the contaminants' velocity vector in Model 3. As shown in Figure9a, the contaminants' flowed from the matrix to the NFs and the fault, because there was a higher permeability, which is consistent with that studied by Zhang[34].Figure9brepresents the concentration curves of the contaminants at the upper boundary.It was observed that the contaminants moved to the upper boundary with a maximum value of 0.2128 mol/m 3 , which was higher than that in Model 1 and Model 2. This shows that the NFs had a significant effect on the amount of contaminants reaching the upper boundary.

Energies 2018 ,Figure 9 .Figure 9 .
Figure 9. (a) Contaminants velocity vector.(b) Contaminants concentration at the upper boundary in Model 3.3.4.Case 4: Alkali Erosion of Rock by High-pH Drilling FluidsSome high-pH (pH = 11-12) drilling fluids are widely used in the drilling process to prevent the wellbore instability, which could cause an alkali erosion of the rock.pH = 12 was chosen in this section to study the influence of alkali erosion on the contaminants' migration.When pH = 12.0, the

4 .
Case 4: Alkali Erosion of Rock by High-pH Drilling Fluids Some high-pH (pH = 11-12) drilling fluids are widely used in the drilling process to prevent the wellbore instability, which could cause an alkali erosion of the rock.pH = 12 was chosen in this section to study the influence of alkali erosion on the contaminants' migration.When pH = 12.0, the initial concentration of OH − was 0.01 mol/L.The relationship between the average reaction rate (k c ) and the pH C OH − was expressed as: k c = 3.03 × 10 −9 C 0.829

Figure 9 .
Figure 9. (a) Contaminants velocity vector.(b) Contaminants concentration at the upper boundary in Model 3.
be obtained from Kang et al.[16].In addition, the increment of rock porosity with time caused by the alkali erosion   c  could be fitted with a third order polynomial, to facilitate the calculation by COMSOL 5.2.

Figure 10 .
Figure 10.The increment of porosity and permeability caused by the alkali erosion of rock.(a) Porosity increment; (b) permeability increment.

Figure 11 .
Figure 11.(a) Concentration distribution of contaminants in Model 4. (b) Contaminants concentration the upper boundary in Model 4.

Figure 11 .
Figure 11.(a) Concentration distribution of contaminants in Model 4. (b) Contaminants concentration at the upper boundary in Model 4.

Figure 11 .
Figure 11.(a) Concentration distribution of contaminants in Model 4. (b) Contaminants concentration at the upper boundary in Model 4.

Figure 13 .
Figure 13.The schematic diagram of geometric model and mesh generation.(a) Geometric model.(b) Mesh generation.Different density grids were adopted in different regions to reduce calculation cost.The NFs, fault and hydraulic fractures regions were described by the refined grids in COMSOL 5.2, and other regions of this model were divided using the standardized grids.

Figure 14 .
Figure 14.Variation in the concentration and accumulation for the three kinds of contaminants.(a) Contaminants concentration at the top of the fault.(b) Contaminants accumulation in aquifers.

Figure 14 .
Figure 14.Variation in the concentration and accumulation for the three kinds of contaminants.(a) Contaminants concentration at the top of the fault.(b) Contaminants accumulation in aquifers.

Figure 14 .
Figure 14.Variation in the concentration and accumulation for the three kinds of contaminants.(a) Contaminants concentration at the top of the fault.(b) Contaminants accumulation in aquifers.

Figure 16 .Figure 17 .
Figure 16.(a) Local distribution of the contaminants velocity vector; (b) variation in velocity at line EF.

Figure 16 .
Figure 16.(a) Local distribution of the contaminants velocity vector; (b) variation in velocity at line EF.

Figure 16 .Figure 17 .
Figure 16.(a) Local distribution of the contaminants velocity vector; (b) variation in velocity at line EF.

Figure 17 .
Figure 17.Variation in the Cd flux and the percentage of the Cd-content, with time.(a) Cd flux; (b) Cd location.

4. 4 .
Analysis of the Influence Factors on the Water Quality of Aquifers 4.4.1.Existence of NFs and Alkali Erosion In this section, three cases (see Table 7) are presented, based on the existence of the NFs and the alkali erosion: (1) Neither NFs nor alkali erosion; (2) with NFs, without alkali erosion; (3) both NFs and alkali erosion.

4. 4 .
Analysis the Influence Factors on the Water Quality of Aquifers 4.4.1.Existence of NFs and Alkali Erosion

Figure 18 .
Figure 18.Variation in the concentration and accumulation for different cases.(a) Contaminants concentration at the top of the fault; (b) contaminants accumulation in aquifers.

Figure 18 .
Figure 18.Variation in the concentration and accumulation for different cases.(a) Contaminants concentration at the top of the fault; (b) contaminants accumulation in aquifers.4.4.2.The Effect of the Number of NFs on the Water Quality of Aquifers

Figure 19 .
Figure 19.Variation in the concentration and accumulation for different numbers of NFs.(a) Contaminants' concentration at the top of the fault; (b) contaminants accumulation in the aquifers.

Figure 19 .
Figure 19.Variation in the concentration and accumulation for different numbers of NFs.(a) Contaminants' concentration at the top of the fault; (b) contaminants accumulation in the aquifers.

Figure 20 .
Figure 20.The increment of porosity and permeability for the different pH values.(a) Porosity increment; (b) permeability increment.

FigureFigure 21 .
Figure 21a,b show the variation curves of the concentration and accumulation of Cd, with time, for the different pH values, respectively.As shown in Figure 21a, the time required for the contaminants migration to the aquifers was 276 years, 260 years, 231 years, 149 years, and 51 years for pH 11.0, 11.5, 12.0, 12.5, and 13.0, respectively.It turned out that higher the pH values, the less time it took for the contaminants to move to the aquifers.In addition, the maximum accumulation of Cd, in the aquifers, was 91.53 mol, 96.76 mol, 106.48 mol, 145.47 mol and 216.48 mol for pH 11.0, 11.5, 12.0, 12.5 and 13.0, respectively.It suggested that the accumulation of heavy metal contaminants in the aquifers increased with an increase of the pH values.

Figure 20 .
Figure 20.The increment of porosity and permeability for the different pH values.(a) Porosity increment; (b) permeability increment.

Figure
Figure 21a,b show the variation curves of the concentration and accumulation of Cd, with time, for the different pH values, respectively.As shown in Figure 21a, the time required for the contaminants migration to the aquifers was 276 years, 260 years, 231 years, 149 years, and 51 years for pH 11.0, 11.5, 12.0, 12.5, and 13.0, respectively.It turned out that higher the pH values, the less time it took for the contaminants to move to the aquifers.In addition, the maximum accumulation of Cd, in the aquifers, was 91.53 mol, 96.76 mol, 106.48 mol, 145.47 mol and 216.48 mol for pH 11.0, 11.5, 12.0, 12.5 and 13.0, respectively.It suggested that the accumulation of heavy metal contaminants in the aquifers increased with an increase of the pH values.

Figure 20 .
Figure 20.The increment of porosity and permeability for the different pH values.(a) Porosity increment; (b) permeability increment.

Figure
Figure 21a,b show the variation curves of the concentration and accumulation of Cd, with time, for the different pH values, respectively.As shown in Figure 21a, the time required for the contaminants migration to the aquifers was 276 years, 260 years, 231 years, 149 years, and 51 years for pH 11.0, 11.5, 12.0, 12.5, and 13.0, respectively.It turned out that higher the pH values, the less time it took for the contaminants to move to the aquifers.In addition, the maximum accumulation of Cd, in the aquifers, was 91.53 mol, 96.76 mol, 106.48 mol, 145.47 mol and 216.48 mol for pH 11.0, 11.5, 12.0, 12.5 and 13.0, respectively.It suggested that the accumulation of heavy metal contaminants in the aquifers increased with an increase of the pH values.

Figure 21 .
Figure 21.Variation in the concentration and accumulation for the different pH values.(a) Contaminants concentration at the top of the fault; (b) contaminants accumulation in the aquifers.

Figure 21 .
Figure 21.Variation in the concentration and accumulation for the different pH values.(a) Contaminants concentration at the top of the fault; (b) contaminants accumulation in the aquifers.

Table 1 .
Characteristics of contaminants from shale gas extraction.

Table 2 .
Summary of important aspects of previous modeling studies.

Table 2 .
Summary of important aspects of previous modeling studies.

Table 3 .
Computational parameters used in model verification.

Table 4 .
Summary of three heavy metal contaminants from the Cambrian Niutitang Formation flow-back liquid.

Table 4 .
Summary of three heavy metal contaminants from the Cambrian Niutitang Formation flow-back liquid.

Table 5 .
Properties of each stratigraphic unit used in model application.

Table 6 .
Computational parameters used in model application.

Table 7 .
Existence of NFs and alkali erosion in different cases.

Table 7 .
Existence of NFs and alkali erosion in different cases.

Table 8 .
The correlation coefficients of the polynomial fit of porosity for the different pH values.

Table 8 .
The correlation coefficients of the polynomial fit of porosity for the different pH values.