Xiaoxu Dong,Wenjing Li,Qun Liu,Zheng Zeng
Research on the Convection-Reaction-Diffusion Model of Contaminants in Fracturing Flowback Fluid in Non-Equidistant Artificial Fractures with Arbitrary Inclination in (3+1)-Dimensional Space-Time
This paper extends the improved fracture equivalent method and then establishes (3+1)-dimensional convection-reaction-diffusion models of contaminants in fracturing flowback fluid in the i-th artificial fracture with arbitrary inclination by comprehensively analyzing the convection effect of fracturing flowback fluid during the flowback process, the diffusion effect of pollutants in fracturing flowback fluid, and the possible chemical reactions between the fracturing fluid and shale matrix. Next, we use a series of transformations and solving methods to solve the established model and obtain semi-analytical solutions of the (3+1)-dimensional convection-reaction-diffusion models. Finally, this paper takes chloride ions as an example to study the concentration changes of pollutants in fracturing flowback fluid in three-dimensional artificial fractures with arbitrary inclination and analyzes the influence of various main control factors on the chloride ion concentration at the inlet end of the i-th artificial fracture with arbitrary inclination.
Research on the Convection-Reaction-Diffusion Model of Contaminants in Fracturing Flowback Fluid in Non-Equidistant Artificial Fractures with Arbitrary Inclination in (3+1)-Dimensional Space-Time

This paper extends the improved fracture equivalent method and then establishes (3+1)-dimensional convection-reaction-diffusion models of contaminants in fracturing flowback fluid in the i-th artificial fracture with arbitrary inclination by comprehensively analyzing the convection effect of fracturing flowback fluid during the flowback process, the diffusion effect of pollutants in fracturing flowback fluid, and the possible chemical reactions between the fracturing fluid and shale matrix. Next, we use a series of transformations and solving methods to solve the established model and obtain semi-analytical solutions of the (3+1)-dimensional convection-reaction-diffusion models. Finally, this paper takes chloride ions as an example to study the concentration changes of pollutants in fracturing flowback fluid in three-dimensional artificial fractures with arbitrary inclination and analyzes the influence of various main control factors on the chloride ion concentration at the inlet end of the i-th artificial fracture with arbitrary inclination.

As a new type of alternative energy, shale gas has a series of advantages such as clean use, wide distribution, low abundance, high production, and long mining life. Therefore, vigorously developing shale gas has become a common choice for countries with rich shale gas resources in the world. Shale gas is mainly bound in tight shale with few pores and fractures. Therefore, large-scale hydraulic fracturing must be carried out for the exploitation of shale gas. The fracturing flowback fluid produced by hydraulic fracturing technology not only contain pressure the highly toxic substances carried by the fracturing fluid itself also contain new pollution elements produced by the chemical reaction between the fracturing fluid and the shale matrix. When the fracturing flowback fluid carries a large amount of highly toxic substances is flowed back to the surface, it will cause surface water pollution and a series of environmental problems if not treated properly. At present, research on flowback of fracturing flowback fluid after hydraulic fracturing in shale gas reservoirs mainly focuses on the methods and models of fracturing flowback fluid flowback and the solution of the models. The following is about shale gas reservoirs and some related investigations on convection, reaction, and diffusion processes of contaminants in fracturing flowback fluid. In 1954, Taylor proposed that in order to reduce the related environmental pollution by formulating relevant decontamination measures in the early stage, it was necessary to study the concentration of contaminants in the flowback fluid during the flowback process of the fracturing flowback fluid. In 1988, Robinson believed that the fracture closure pressure was greatly affected by the size of the oil nozzle. On the basis of analyzing the fracture closure pressure under the small nozzle size, a theoretical model of fracturing flowback fluid with small displacement was established. Its main feature was to control the speed of the flowback. In 2004, according to the seepage mechanics and the superposition of potentials, Xu comprehensively considered factors such as the number of fractures, the spacing of fractures, the length of fractures, the symmetry of fracture distributions, and the azimuth angle of fractures. A new model for calculating the productivity of horizontal wells after fracturing is established, and the analytical solution of the new model is given. In 2006, Luo gave the equations for calculating the flow velocity in the fractures generated by hydraulic fracturing in one-dimensional and two-dimensional cases, and at the same time, without considering the influence of the compressibility of the fracturing fluid in the fractures on the flow rate. According to the equilibrium volume principle, the continuity equation of fluid flow in fractures is given. In 2007, Jiang considered the diffusion, adsorption, and desorption of organic pollutants and microbial degradation in the saturated zone and established a mathematical model for the affinity of groundwater flow and pollutant migration and studied the influence of parameters on the migration law. In 2007, Schmidt solved the exact solution of the two-component equations of the reaction–diffusion equations with the method of linearly deterministic equations. The limit differential equations of the second-order evolution equation system are constructed. In 2009, Moitsheki et al. studied the transient coupling of partial differential equations describing the fluid flow and mass transfer of pollutants in a cylindrical coordinate system. The symmetry analysis of the coupled partial differential equation system is carried out, and some special Lie algebras are obtained for any special choice of constants and source terms. For all admissible symmetries, the optimal system is formed. In 2012, Xu selected three representative pollutants to conduct a numerical simulation study on the migration and pollution remediation laws of organic matter and heavy metals in groundwater. And the adsorption and desorption processes, chemical reactions, and biodegradation functions of pollutants in porous media are added as the source or sink term of the convection–dispersion equation into the convection–dispersion equation of pollutant migration in porous media. In 2014, Abbasi et al. estimated the fracture permeability and fracture volume using a simple straight-line technique. Most of the results show that these two parameters show a high linear correlation on the double logarithmic curve. In 2015, Liu established a mathematical model for the migration of pollutants in porous media based on the kinetics of pollutant migration in porous media. This study revealed the intrinsic relationship of convection, dispersion, adsorption, and other behaviors of pollutants in porous media. In 2016, Bernardi et al. described organic pollutants in mixed reaction–diffusion systems, proposed a mixed variational form, and reviewed its good possibilities. Then, considering the spectral discretization of this problem, an almost optimal error estimate is established. In 2016, Bhadane et al. adopted a new integral transformation known as the Elzaki transformation to solve a model consisting of two convection–diffusion equations. In 2016, Park established a two-dimensional hydrodynamic model and a pollutant mixing model and considered flow and migration phenomena. The governing equations of the pollution transport model provide a more accurate model for the prediction of dispersion. In 2017, Li et al. gave an analytical solution to the convection–diffusion equations through the two-dimensional Fourier transformation and the inverse Fourier transformation. In 2019, Liu et al. established a mathematical model for hydraulic fracturing flowback fluid in shale gas wells, which provided useful information for the complexity of the fracture network after fracturing according to the dynamic changes of the flowback fluid. In 2019, Prakash discussed a complete classification of invariant subspaces for a class of generalized nonlinear convection-reaction-diffusion equations, which are widely used to describe many natural phenomena, such as fluid convection in porous media and chemical reaction–diffusion equations. In 2022, Dong et al. proposed the improved fracture equivalent method and solved the convection-reaction-diffusion models of contaminants in fracturing flowback fluid in two-dimensional non-equidistant artificial fractures with arbitrary inclination. At present, the improved fracture equivalent method is only applicable to the case where the fractures are distributed in the two-dimensional plane rectangular coordinate system, and the numerical solution of the three-dimensional convection-reaction-diffusion model is mostly obtained by using numerical methods. Therefore, in order to study the convection, reaction, and diffusion processes in three-dimensional non-equidistant artificial fractures with arbitrary inclination, the main tasks of this paper are as follows. (1) The improved fracture equivalent method is extended to three-dimensional space rectangular coordinate systems, which is originally only applicable to two-dimensional plane rectangular coordinate systems. The improved fracture equivalent method after expansion is effectively applicable to artificial fractures with arbitrary inclination in three-dimensional space. (2) This paper establishes and solves convection-reaction-diffusion models of contaminants in fracturing flowback fluid in three-dimensional non-equidistant artificial fractures with arbitrary inclination. A series of transformations and scientific methods are used to obtain semi-analytical solutions for (3+1)-dimensional convection-reaction-diffusion models.

In previous research, we have discussed in detail the convection, reaction, and diffusion mechanisms of contaminants in fracturing flowback fluid in fractures. Based on the convection, reaction, and diffusion mechanisms, this section establishes (3+1)-dimensional convection-reaction-diffusion models of contaminants in fracturing flowback fluid in i-th artificial fracture with arbitrary inclination. The convection–diffusion equation can be obtained asThe migration of contaminants also includes adsorption and desorption processes, chemical reaction processes, and biodegradation processes. These processes can be added to eq 1 as the source or sink term of the convection–diffusion equation, which isThe first-order chemical reaction kinetics equation (reaction rate proportional to pollutant concentration) is widely used in the problem of pollutant transport to describe the variation of pollutant concentration over time during chemical and biological reactions. The expression of the first-order chemical reaction kinetics equation iswhere ζ is the first-order reaction rate constant and C is the concentration of contaminants in the liquid phase. When ζ < 0, it indicates a continuous decrease in the concentration of contaminants in the liquid phase, but the opposite is true. Initial condition: Boundary conditions: Dirichlet boundary condition: Neumann boundary condition: Cauchy boundary condition: Taking the i-th artificial fracture as the research object, the convection-reaction-diffusion schematic diagrams of the contaminants in the fracturing flowback fluid on the left and right wings of the artificial fractures are shown in Figures 1 and 2. Based on the study of convection-reaction-diffusion model of fluid, the (3+1)-dimensional convection-reaction-diffusion models of contaminants in fracturing flowback fluid under three outer boundary conditions in the right wing of the i-th non-equidistant artificial fracture with arbitrary inclination are as follows:In the same way, the (3+1)-dimensional convection-reaction-diffusion models of contaminants in fracturing flowback fluid under three outer boundary conditions in the left wing of the i-th non-equidistant artificial fracture with arbitrary inclination are as follows:

In order to obtain the tip coordinates (xfil, yfil, zfil) and (xfir, yfir, zfir) of non-equidistant artificial fractures with arbitrary inclination in three-dimensional space, this paper further studies the improved fracture equivalent method considering the asymmetry of the left and right wings of the fractures. The improved fracture equivalent method proposed by Dong is only applicable to the case where the left and right fractures are distributed in the two-dimensional Cartesian coordinate system. However, the actual distributions of fractures in three-dimensional space are more complicated. In order to better describe all possible distribution of artificial fractures in three-dimensional space, this section further extends the improved fracture equivalent method on the basis of the original foundation. The extended improved fracture equivalent method is more suitable for the actual situation. Assume that the direction of the horizontal wellbore is the positive direction of the x-axis, and the direction perpendicular to the horizontal wellbore and pointing upward is the positive direction of the z-axis. The angle between the left wing of the fracture and the positive direction of the x-axis is φfil(i)(0 ≤ φfil(i) ≤ 180°), and the angle between the left wing of the fracture and the positive direction of the z-axis is θfil(i)(0 ≤ θfil(i) ≤ 180°). The angle between the right wing of the fracture and the positive direction of the x-axis is φfir(i)(0 ≤ φfir(i) ≤ 180°), and the angle between the right wing of the fracture and the positive direction of the z-axis is θfir(i)(0 ≤ θfir(i) ≤ 180°). Consistent with the improved fracture equivalent method, the left and right wings of the non-equidistant artificial fracture with arbitrary inclination are divided into n equal parts, and each part is regarded as a point sink to study, and the coordinates of each point sink are the midpoint coordinates of each equal part. In the spatial Cartesian coordinate system, the schematic diagrams of the azimuth and fracture bisector of the horizontal wellbore of artificial fractures are shown in Figures 3 and 4. The 16 schematic diagrams in Figures 3 and 4 show the distributions of all possible artificial fractures in three-dimensional space. Let the equivalent length of the j-th point sink on the right (left) wing of the i-th artificial fracture with arbitrary inclination on the perforation section of the horizontal wellbore converge to the interface between the horizontal wellbore and the right (left) wing of the i-th artificial fracture with arbitrary inclinationUnder the 16 different fracture distributions, the coordinates of the j-th point sink on the right wing of the i-th artificial fracture with arbitrary inclination on the perforation section of the horizontal wellbore areUnder the 16 different fracture distributions, the coordinates of the j-th point sink on the left wing of the i-th artificial fracture with arbitrary inclination on the perforation section of the horizontal wellbore are:From eqs 15 and 16, the coordinates of the left and right wing tips of the i-th artificial fracture with arbitrary inclination can be obtained as follows:where 0 < α(i), β(i) < π, and i = 1,2, ···, m.

In this section, we give the solutions of the convection-reaction-diffusion models established in Section 2 for the contaminants in the fracturing flowback fluid in artificial fractures with arbitrary inclination.

The following coordinate transformation is performed for the definite solution problems (8–10) of the partial differential equationThen, the (3+1)-dimensional definite solution problems of the partial differential equation are transformed into the following (1+1)-dimensional definite solution problems of the partial differential equation.where Λ(fir) = [sinθ(fir) cosφ(fir) + sinθ(fir) sinφ(fir) + cosθ(fir)] and the direction of the radial direction r is the direction of the centerline of the fracture. Next, introducing dimensionless variablesDimensionless parameters list:According to eq 23, it can be known thatSubstituting eqs 26–28 into the definite solution problems (20–22) of the partial differential equation, thenFor the definite solution problems (29–31) of the partial differential equation, taking the following Laplace transformations on the dimensionless time variable TThen, eqs 29–31 are transformed into the following definite solution problems of ordinary differential equation.The solutions of definite solution problems (34–36) in Laplace space arewhere C̅DDfir(R,S), C̅NDfir(R,S), and C̅CDfir(R,S) are the Laplace spatial concentrations of contaminants in the right wing of the i-th artificial fracture with arbitrary inclination under the Dirichlet, Neumann, and Cauchy boundary conditions. The Laplace spatial concentrations of contaminants in the fracturing flowback fluid in the right wing of the i-th non-equidistant artificial fracture with arbitrary inclination are Using the modified Gaver–Stehfest numerical inversion method, the concentrations of contaminants in the fracturing flowback fluid in the right wing of the i-th non-equidistant artificial fracture with arbitrary inclination in real space can be obtained aswhere In real space, the concentrations of contaminants in the fracturing flowback fluid at the interface between the right wing of the i-th artificial fracture with arbitrary inclination and the horizontal wellbore are

The following coordinate transformation is performed for the definite solution problems (11–13) of the partial differential equation Then, the (3+1)-dimensional definite solution problems of the partial differential equation are transformed into the following (1+1)-dimensional definite solution problems of the partial differential equation.where Λ(fil) = [sinθ(fil) cosφ(fil) + sinθ(fil) sinφ(fil) + cosθ(fil)] and the direction of the radial direction r is the direction of the centerline of the fracture.Next, introducing dimensionless variablesDimensionless parameters list:According to eq 48, it can be known thatSubstituting eqs 51–53 into the definite solution problems (45–47) of the partial differential equation, thenFor the definite solution problems (54–56) of the partial differential equation, taking the following Laplace transformations on the dimensionless time variable TThen, the definite solution problems (54–56) of the partial differential equation are transformed into the following definite solution problems of ordinary differential equationThe definite solution problems (59–61) of ordinary differential equation are similar to that of definite solution problems (34–36), and the solutions of definite solution problems (59–61) in Laplace space arewhere C̅DDfil(R,S), C̅NDfil(R,S), and C̅CDfil(R,S) are the Laplace spatial concentrations of contaminants in the left wing of the i-th artificial fracture with arbitrary inclination under the boundary Dirichlet, Neumann, and Cauchy conditions. The Laplace spatial concentrations of contaminants in the fracturing flowback fluid in the left wing of the i-th non-equidistant artificial fracture with arbitrary inclination areUsing the modified Gaver–Stehfest numerical inversion method, the concentrations of contaminants in the fracturing flowback fluid in the left wing of the i-th non-equidistant artificial fracture with arbitrary inclination in real space can be obtained aswhere It can be seen that in real space, the concentrations of contaminants in the fracturing flowback fluid at the interface between the left wing of the i-th non-equidistant artificial fracture with arbitrary inclination and the horizontal wellbore are

In this section, the Dirichlet boundary condition is used as an example to calculate the concentration of contaminants at the interface between the right wing of the i-th artificial fracture and the horizontal wellbore, in which φfir(i) is the inclination with the positive direction of the x-axis, and θfir(i) is the inclination with the positive direction of the z-axis.

Letwhere afir, bfir, and hfir are known constants. Because ofthen Equation 69 is transformed intoTaking dimensionless transformation on eq 72, thenSubstituting eq 73 into eq 37, the concentration of contaminants at the interface between the right wing of the i-th artificial fracture and the horizontal wellbore under the Dirichlet boundary condition is In real space, the concentration of contaminants at the interface between the right wing of the i-th artificial fracture and the horizontal wellbore is Taking chloride ions as an example, the relevant parameter values of the eq 74 are shown in Table 1. According to eqs 67, 74, and 75, the corresponding programs are compiled, and the obtained concentration curve is shown in Figure 5. The composition of fracturing flowback fluid is complex, which includes about 10–70% of fracturing fluid, heavy metals, and soluble solids. Besides, the soluble solids mainly include Na+, K+, and chloride, and the total dissolved solids content is extremely high. At the same time, the longer the residence time of the fracturing fluid in fractures, the higher the content of soluble solids in the flowback fluid, which means that the chloride ion concentration in the flowback fluid gradually increases with time during the initial stage of flowback. Artificial fractures are generally considered as large fractures with stronger conductivity than natural fractures. Therefore, when the reaction between fracturing fluid and shale matrix gradually reaches a certain saturation state, the chloride ion concentration rapidly decreases under corresponding fracture parameter settings. The concentration curve presents a single-peak-type curve. Based on the researches of Yu and Xu, it can be seen that the model is correct.

On this basis, we studied the influence of the flowback velocity vfir of the fracturing flowback fluid, the first-order reaction coefficient ζfir of the pollutant concentration, the molecular diffusion coefficient D of contaminants in the flowback fluid, and the length Lfir of the right wing of the i-th artificial fracture on the change of the concentration of chloride ions in fracturing flowback fluid in the right wing of i-th non-equidistant artificial fractures with arbitrary inclination. According to Figure 6, it can be seen that under the Dirichlet outer boundary condition, while other parameters remain unchanged, the dimensionless peak concentration of chloride ions at the inlet of the right wing of the i-th artificial fracture decreases with the linear increase of the flowback velocity vfir of the fracturing fluid in the right wing of the i-th artificial fracture. The reason is that when the length of the i-th artificial fracture remains unchanged, as the flowback velocity vfir of the fracturing fluid increases in the i-th artificial fracture, the time that the fracturing fluid stays in the fracture becomes shorter. Therefore, the reaction time between fracturing fluid and shale matrix becomes shorter, ultimately leading to a decrease in the peak concentration of chloride ions in the flowback fluid at the right wing inlet of the i-th artificial fracture. According to Figure 7, it can be seen that under the Dirichlet outer boundary condition, when other parameters remain unchanged, with the linear increase of molecular diffusion coefficient D, the dimensionless peak concentration of chloride ion at the inlet of the right wing of the i-th artificial fracture increases. At the same time, it can be seen that the change time of chloride ion concentration is earlier, and the growth time of chloride ion concentration in the flowback stage is longer. The reason is that the fracturing fluid will diffuse from the horizontal wellbore into the fracture. As the molecular diffusion coefficient D increases, it means that the diffusion degree of chloride ions in the fracture increases. Therefore, the increment of chloride ions generated by the diffusion effect increases, and it leads to an increasing trend in the peak concentration of chloride ions at the right wing inlet of the i-th artificial fracture. According to Figure 8, it can be seen that under the Dirichlet outer boundary condition, while other parameters remain unchanged, as the first-order chemical reaction coefficient ζfir linearly decreases, the dimensionless peak concentration of chloride ions at the inlet of the right wing of the i-th artificial fracture decreases. The reason is that as the first-order chemical reaction coefficient ζfir in the right wing of the i-th artificial fracture decreases, the amount of chloride ion reduction in the right wing of the i-th artificial fracture is greater due to the chemical reaction between the fracturing fluid and shale matrix. This ultimately leads to a decrease in the dimensionless peak concentration of chloride ions at the inlet of the right wing of the i-th artificial fracture. According to Figure 9, while other parameters remain unchanged, as the length Lfir of the right wing of the i-th artificial fracture linearly increases, the dimensionless peak concentration of chloride ion at the inlet of the right wing of the i-th artificial fracture decreases. The reason is that as the length Lfir of the right wing of the i-th artificial fracture increases, the path of the fracturing fluid flowing back from the right wing of the i-th artificial fracture into the horizontal wellbore becomes longer. This results in a time lag in the beginning of the change in chloride ion concentration. At the same time, the longer the length of the right wing of the i-th artificial fracture, the greater the space occupied by the fracture, and the more fracturing fluid enters the right wing of the i-th artificial fracture. Therefore, the chemical reaction between the fracturing fluid and the shale matrix is more intense, and the first-order chemical reaction coefficient is less than 0. This ultimately led to a decrease in the dimensionless peak concentration of chloride ions at the inlet of the right wing of the i-th artificial fracture.

(1) This paper further extends the improved fracture equivalent method to three-dimensional space and depicts 16 distributions of left and right wings of non-equidistant artificial fractures with arbitrary inclination in three-dimensional space. The improved fracture equivalent method after expansion is more suitable for the actual distribution of artificial fractures generated by hydraulic fracturing in a shale matrix. (2) In this paper, a series of transformations and scientific methods, such as coordinate transformation, dimensionless transformation, Laplace transformation, undetermined coefficient method, Gaver–Stehfest numerical inversion method, and special functions, are used to solve the convection-reaction-diffusion models of pollutants in flowback fluid in three-dimensional non-equidistant artificial fractures with arbitrary inclination. The semi-analytical solutions of the pollutant concentration at the interface between the left and right wings of the i-th non-equidistant artificial fracture with arbitrary inclination and the horizontal wellbore are obtained.