Mathematical Modeling of Hydraulic Fracture Formation and Cleaning Processes

: The effectiveness of the hydraulic fracturing procedure is crucially dependent on the stage of fracture planning and design. Forecasting fracture behavior in rock formations characterized by non-uniform toughness is a serious challenge. In the present paper, a planar-3D model considering the rock’s non-uniform fracture toughness has been developed for the uneven propagation of a hydraulic fracture. The series of numerical experiments were designed to study the effect of inhomogenous fracture toughness. The results show that the fracture toughness contract signiﬁcantly controls the overall direction of fracture propagation, and a combination of toughness contrast and the proportion between the pay zone and barrier zone determine the fracture proﬁle: from almost circular with or without a pair of narrow wedges when the proportion is small to almost rectangular otherwise. This paper also discusses the process of cleaning a fracture from hydraulic fracturing ﬂuid by oil. Using numerical modeling on the basis of the constructed mathematical model, a relationship is established between the quality of hydraulic fracture cleaning and the geometrical parameters of the fracture and the region ﬁlled with the hydraulic fracturing ﬂuid. The results of numerical experiments show that while fracturing ﬂuid is more viscous than oil, the length of the fracture has a greater inﬂuence on the cleaning process than the viscosity of the fracturing ﬂuid.


Introduction
With the increasing demand for fossil fuels (oil and gas) and the depletion of conventional oil and gas reservoirs, unconventional reservoirs have attracted great attention and now occupy a large proportion of the hydrocarbon supply [1,2]. Previously, unconventional reservoirs were difficult for exploitation due to their extremely low permeability and conductivity [3,4]. Hydraulic fracturing is an efficient technique for improving unconventional reservoirs' permeability and conductivity via hydraulic fractures [5]. In the hydraulic fracturing design fracture geometry (length, width and height) is a critical target, which is determined by two "controllable" variables: hydraulic fluid viscosity and wellbore pumping rate, and four "natural" variables: in situ stress, reservoir's stiffness (Young's modulus and Poisson's ratio), fluid loss (porosity and permeability) and fracture toughness [6]. Hydraulic fluid viscosity and wellbore pumping rate are controlled artificially [7][8][9], so they can be easily studied and their working mechanism is relatively clear. The remaining four "natural" variables are mainly controlled by geological conditions, which have strong randomness and complexity. Therefore, the attention of researchers was mainly focused on studying the effects of these four "natural" variables on fracture geometry and propagation. The sustainability of the hydraulic fracturing procedure is crucially dependent on the stage of fracture planning and design. Forecasting fracture behavior in rock formations characterized by non-uniform toughness is a serious challenge.
Sedimentary rocks around the reservoir usually exhibit layered homogeneity due to the deposition and compaction processes [10,11]. Different layers have their own stiffness, porosity and permeability, fracture toughness as well as in situ stress. This situation causes significant contrasts of the four "natural" variables between layers, especially between the pay zone and the adjacent layers. In the early stages of fracturing studies, Warpinski [12,13], conducted experiments to determine which of the two parameters (stiffness contrast and in situ stress contrast) is predominant for controlling hydraulic fracture containment. Experiments demonstrated that the in situ stress contrast is the most important factor controlling fracture height, while the stiffness contrast has little effect. This conclusion is consistent with numerical and semi-analytical results in the paper [14]. In contrast to conclusions in [12][13][14], papers [6,15] gave quantitative descriptions of the effects of formation stiffness (Young's modulus and Poisson's ratio) contrast on hydraulic fracture propagation. In [16,17], pseudo-3D models were presented for hydraulic fracture growth in a layered formation with contrasts in both stiffness and in situ stress. In addition, some attention was paid to studying the effects of fluid loss (porosity and permeability) on hydraulic fracture propagation. Warpinski et al. [13] presented a qualitative conclusion that permeability and pore pressure have an important effect on hydraulic fracture containment. Hanson et al. [18] developed theoretical and numerical models and indicated that the increase in the pore pressure tends to reduce the tendency for crack extension, and the fracture growth could be impeded in areas where the permeability is large. In papers [19,20], the models for hydraulic fracture propagation in an inhomogeneous poroelastic medium were proposed for demonstrating the influence of reservoir permeability on the dynamics of fracture propagation. Gao et al. [21] established a 2D fluid-solid coupled model to study the hydraulic fracture propagation in a non-uniform pore pressure field. The results indicated that the distributions of in situ stress and pore pressure are gradually altered with the fracture propagation, and these alterations inversely affect subsequent fracture growth. Furthermore, different models and algorithms were proposed for the combined effects of the four "natural" variables on hydraulic fracture propagation in a layered medium [17,[22][23][24][25].
In the above descriptions of the effects of the four "natural" variables (in situ stress, stiffness, fluid-loss and fracture toughness) on hydraulic fracture propagation in a multilayered medium there is a common feature: ignoring or weakening the influence of fracture toughness parameter, because most researchers think that fracture toughness has a negligible effect on hydraulic fracturing, and it shows a relatively small variation range in unconventional reservoirs [26]. In fact, relatively accurate distributions of fracture toughness in the layered media of reservoirs are very rare, because the measurement methods for fracture toughness in conditions of fracture growing under complex geological conditions are still immature [27][28][29][30]. Thiercelin et al. [31,32] performed experiments and numerical studies on the influence of fracture toughness in a homogeneous reservoir as well as on the influence of fracture toughness contrast in a layered medium. The results show that fracture toughness has a significant effect on fracture propagation not only in the homogeneous formation, but also in the layered medium. When the contrast in fracture toughness is sufficiently high, the arrest or diversion of fracture may be obtained. It was also pointed out in [32], that fracture toughness is a critical mechanical property influencing hydraulic fracture propagation, especially in the cases where the in situ stress contrast is small, the fluid is of low viscosity, and the fracture is relatively small.
Since fracture toughness is not the main factor in the above four "natural" variables, and the data of fracture toughness in the layered medium are difficult to measure, there are very few studies of the effect of fracture toughness contrast. In order to better understand the mechanism of the effect of non-uniform fracture toughness on hydraulic fracturing, especially in the early stage when the fracture is small, in this paper we established a plane-3D model considering an inhomogeneous fracture toughness for hydraulic fracture Energies 2022, 15, 1967 3 of 35 propagation. In this model, the three-dimensional deformation of the fracture is coupled with the two-dimensional flow of hydraulic fluid, and the two-directional propagation in the fracture plane is described in the framework of linear elastic fracture mechanics in the cylindrical coordinate system. Besides, a very small-time step is applied for capturing the uneven propagation on different locations of the fracture front. Furthermore, many numerical experiments have been carried out for quantitative analysis of the effect of fracture toughness on hydraulic fracture propagation.

Mathematical Model for Hydraulically Driven Fracture Propagation
Consider a vertical hydraulic fracture across a horizontal wellbore in a layered medium, as illustrated in Figure 1. For clarity of discussion on the effect of inhomogeneity of fracture toughness, the medium is assumed here to be impermeable and other properties of the medium, Young's modulus E, Poisson's ratio ν and in situ stress σ 0 are assumed uniform and homogeneous, except the fracture toughness. Hydraulic fluid is injected at a constant rate Q 0 and is assumed to follow a Newtonian rheology with a viscosity µ. a plane-3D model considering an inhomogeneous fracture toughness for hydraulic fracture propagation. In this model, the three-dimensional deformation of the fracture is coupled with the two-dimensional flow of hydraulic fluid, and the two-directional propagation in the fracture plane is described in the framework of linear elastic fracture mechanics in the cylindrical coordinate system. Besides, a very small-time step is applied for capturing the uneven propagation on different locations of the fracture front. Furthermore, many numerical experiments have been carried out for quantitative analysis of the effect of fracture toughness on hydraulic fracture propagation.

Mathematical Model for Hydraulically Driven Fracture Propagation
Consider a vertical hydraulic fracture across a horizontal wellbore in a layered medium, as illustrated in Figure 1. For clarity of discussion on the effect of inhomogeneity of fracture toughness, the medium is assumed here to be impermeable and other properties of the medium, Young's modulus , Poisson's ratio and in situ stress are assumed uniform and homogeneous, except the fracture toughness. Hydraulic fluid is injected at a constant rate and is assumed to follow a Newtonian rheology with a viscosity . Figure 1. A planar-3D hydraulic fracture across a horizontal wellbore in a medium with inhomogeneous fracture toughness.

Governing Equations
For establishing the planar-3D model for hydraulic fracture propagation in an impermeable medium we consider: the three-dimensional deformation of fracture, the two-dimensional flow of hydraulic fluid in the fracture, and the two directional propagations in the inhomogeneous medium. This corresponds to the planar 3D approach described in the introductory section, in which the three-dimensional problem of hydraulically driven fracture propagation in a non-uniform medium is treated under the assumption that the fracture propagates within a definite plane and is symmetrical in respect to that plane.

Governing Equations
For establishing the planar-3D model for hydraulic fracture propagation in an impermeable medium we consider: the three-dimensional deformation of fracture, the twodimensional flow of hydraulic fluid in the fracture, and the two directional propagations in the inhomogeneous medium. This corresponds to the planar 3D approach described in the introductory section, in which the three-dimensional problem of hydraulically driven fracture propagation in a non-uniform medium is treated under the assumption that the fracture propagates within a definite plane and is symmetrical in respect to that plane. This is an assumption, which is supported by the majority of test cases, nevertheless, it is not always the case. This assumption provides the possibility for a serious simplification of the mathematical model regarding the fracture propagation and fluid flow within the frame of a 2-D approximation thus reducing the dimension of the problem and reducing essentially the computer simulation time. On the other hand, the stress-strained state of the medium surrounding the fracture is treated in a full scale 3-D statement. Given the elastic properties, toughness and confining stress of the medium, the fluid properties, the injection rate and the initial fracture, the solution of fracture propagation is composed of the fluid pressure p f (x, y, t), the fracture width, namely fracture opening w(x, y, t) and the projection area Ω(t) of the fracture on the plane xoy.

Elastic Deformation
For the planar fracture propagation, namely mode I fracture, the fluid pressure in the fracture can be expressed by the following integral equation [33] p(x, where p(x, y, t) is the net pressure in the fracture, E = E/ 1 − ν 2 is the plane strain elastic modulus.

Fluid Flow
According to the lubrication theory and Poiseuille's law, the generalized continuity equation without considering the fluid loss is where → u is the fluid velocity inside the fracture, ∇p is the net pressure gradient.

Boundary Conditions
Since the fluid is injected into the fracture with a constant rate Q 0 , the fluid flux at the wellbore is formulated as where w 0 and u 0 are the fracture width and fluid velocity at the wellbore, respectively, and r 0 is the wellbore radius. According to the linear elastic fracture mechanics and the assumption of zero fluid lag, the fracture width and the fluid flux at the fracture front equal zero, namely where s c is the boundary curve of the projection area Ω(t) of the fracture on the plane xoy. The additional boundary condition at the fracture front is given by the typical fracture propagation criterion: where K I (s c , t) is the stress intensity factor with respect to parameters time and fracture front coordinates, K IC (s c ) is the fracture toughness with respect to fracture front coordinates because of inhomogeneity.

Initial Conditions
It is generally known that fracture initiation is a complex process [34,35], which is not the focus of our model. To simplify the model, we make an assumption that at the initial time of simulation (t 0 = 0) the fracture has a small penny-like shape with a radius R 0 and is opened by fluid pressure with a uniform distribution p 0 (Figure 2), therefore, the cylindrical coordinate system rθz is chosen for better describing the penny-shape fracture Energies 2022, 15, 1967 5 of 35 evolution, especially the two-dimensional fluid flow. In the cylindrical coordinate system, for formulating the layered homogeneous fracture toughness of formation we assume that fracture toughness is only related to the angular coordinate θ, namely and K IC (θ) = K IC (θ + π) (0 ≤ θ ≤ π). For example, in Figure 2 the thick red line drawn across the point source represents a homogeneous rock layer with a uniform fracture toughness: K IC (0) = K IC (π).
Energies 2022, 15, x FOR PEER REVIEW 5 of 36 time of simulation ( = 0) the fracture has a small penny-like shape with a radius and is opened by fluid pressure with a uniform distribution ( Figure 2), therefore, the cylindrical coordinate system is chosen for better describing the penny-shape fracture evolution, especially the two-dimensional fluid flow. In the cylindrical coordinate system, for formulating the layered homogeneous fracture toughness of formation we assume that fracture toughness is only related to the angular coordinate , namely and ( )= ( + ) (0 ≤ ≤ ). For example, in Figure 2 the thick red line drawn across the point source represents a homogeneous rock layer with a uniform fracture toughness: (0)= ( ).

Model in Cylindrical Coordinate System
The governing Equation (1) of fracture deformation can be easily transformed into the cylindrical coordinate system with the help of the relationship between two systems: = cos , = sin . It should be noted that we can solve this equation in the Cartesian coordinate system first, and then transform the final results into the cylindrical system instead of directly substituting the transformation relationship into Equation (1) for solving it in the system . The purpose of this is to simplify the calculation greatly, because Equation (1) is a special hyper singular integral expression, which can be calculated easier in the Cartesian system than in the polar system . The governing Equation (2) can be transformed into the polar system as

Model in Cylindrical Coordinate System
The governing Equation (1) of fracture deformation can be easily transformed into the cylindrical coordinate system with the help of the relationship between two systems: x = r cos θ, y = r sin θ. It should be noted that we can solve this equation in the Cartesian coordinate system xoy first, and then transform the final results into the cylindrical system instead of directly substituting the transformation relationship into Equation (1) for solving it in the system roθ. The purpose of this is to simplify the calculation greatly, because Equation (1) is a special hyper singular integral expression, which can be calculated easier in the Cartesian system xoy than in the polar system roθ. The governing Equation (2) can be transformed into the polar system roθ as where µ = 12µ. In this equation the radial and angular fluid flows in the fracture are considered: The boundary conditions at the wellbore and at the fracture front are rewritten as where R(θ) is the fracture radius with respect to θ, because the fracture has different radii at different radial directions.

Numerical Scheme
The closed equation system of this model is composed of Equations (1), (7), (9) and (10). The unknowns net pressure p(r, θ, t) and fracture opening w(r, θ, t) are coupled by governing Equations (1) and (7) in so complex form that only the implicit scheme can be used to solve this system. Furthermore, it should be noted that Neumann boundary conditions (Equations (9) and (10)) require a special consideration [36]: the boundary pressure at the fracture front p(R, θ, t) as an estimated known variable has to enter into the system for a solution, and then the accurate numerical solution of p(R, θ, t) will be obtained by iteration according to the law of conservation of fluid volume:

Fracture Discretization
The fracture plane is discretized with a fixed polar grid using triangular and trapezoidal elements ( Figure 3). ∆r and ∆θ are the grid step lengths in the r and θ directions respectively. The main physical parameters (location coordinates, fracture opening and fluid pressure) are evaluated at the element centers of gravity. Besides, the elements with the same angular coordinate θ j have the same fracture toughness K j IC , just as the elements passed by the red dotted line in Figure 3. The boundary elements marked with black points contain an extra parameter, namely the stress intensity factor K j I . If at a certain moment there are some boundary elements satisfying the condition K j I ≥ K j IC , the fracture front where these boundary elements are located would propagate one element forward, see the new elements with red dotted curves at directions θ j+1 and π + θ j+1 .

Discretization of Equation System
As already mentioned in Section 2.2, Equation (1) is solved in the Cartesian system first, and then the obtained results are transformed into the polar system. This can greatly simplify the calculation process, so Equation (1) is discretized as where ( , ) are coordinates of any fracture element in Cartesian system, ( , ) are coordinates of each fracture element in the Cartesian system. is the element number in the direction of the farthest fracture element away from the wellbore, = 2 ⁄ . The stiffness coefficient ( , , , ) is formulated as  . The discretized polar grid for the fracture. Black points-boundary grid nodes of the fracture front, red points-internal grid nodes, blue points-nodes of gravity center of elements and the green point is the wellbore. r i , θ j and r i , θ j+1 are the elements' coordinates. Subscripts "i " and "j " are the element number in the r and θ directions, in which j = 1, 2, 3, · · · , n, and n = 2π/∆θ.

Discretization of Equation System
As already mentioned in Section 2.2, Equation (1) is solved in the Cartesian system first, and then the obtained results are transformed into the polar system. This can greatly simplify the calculation process, so Equation (1) is discretized as where (x, y) are coordinates of any fracture element in Cartesian system, x i , y j are coordinates of each fracture element in the Cartesian system. m is the element number in the r direction of the farthest fracture element away from the wellbore, n = 2π/∆θ. The stiffness coefficient I x, y, x i , y j is formulated as where A x i , y j is the area of fracture element x i , y j . The coordinates of fracture elements in the two kinds of coordinate systems can be transformed into each other with the help of the following relationship: According to the finite volume method, Equation (7) for fluid flow can be discretized as follows where ∆w = w i,j,k+1 − w i,j,k represents the increment of fracture opening of the element r i , θ j at the moment t = t k+1 . Subscript "k" is the time step number. Boundary conditions (9) and (10) are transformed as where subscripts "1 − 1/2" and "m + 1/2" in the r direction represent the inner boundary node at the wellbore and the outer boundary node on the fracture front, respectively. For the utilization of Equation (17), the pressure at the fracture front p m+1/2,j should enter into Equations (12), (15)- (17) as an estimated known parameter [36], and then be solved by iteration according to Equation (11). The stress intensity factor K I in the three-dimensional model can be solved by its definition: the limit of the product of the stress near the fracture tip and the root of the distance from the fracture tip. Bui [37], proposed an equivalent method for solving K I in the three-dimensional model with a precondition that the discontinuous displacement near the fracture tip is known. Its discrete form is where w m,j and r m are the fracture opening and radial coordinate of the nearest element away from the fracture front in each direction θ = θ j respectively. R j is the fracture radius in the direction θ = θ j .

Implicit Algorithm
The discretized model consists of Equations (12), (15)- (17) with the tip condition of stress intensity factor (18). Firstly, substituting Equation (12) into (15)-(17), we obtain an implicit nonlinear equation system for the unknowns w i,j,k+1 with the highest fourth power, which can be solved by the Newtonian iteration method. Then, the obtained values of w i,j,k+1 must satisfy the law of fluid volume conservation (11): if condition (11) is satisfied, proceed into the second step; otherwise, further iteration is needed for w i,j,k+1 until the condition of fluid volume conservation (11) is satisfied.
Secondly, the stress intensity factors K I θ j in all directions θ can be calculated (Equation (18)) using the obtained values of w i,j,k+1 . If all stress intensity factors satisfy the condition K I θ j < K IC θ j , it means that the values of w i,j,k+1 are the final solutions of the model at the moment t = t k+1 and the third step for solving p i,j,k+1 starts, otherwise, it means that fracture front where stress intensity factors satisfy K I θ j ≥ K IC θ j propagates, and the fracture planar area with the fracture elements is updated. Therefore, the entire Energies 2022, 15, 1967 9 of 35 calculation area for the model is updated, and the calculation needs to return to the first step for new values of w i,j,k+1 , until all stress intensity factors K I θ j satisfy the condition (18). The final values of w i,j,k+1 as well as the new fracture front satisfying the conditions (11) and (18) are the numerical solutions of the model at the moment t = t k+1 .
Thirdly, the obtained numerical values of w i,j,k+1 are inversely substituted into Equation (12) for the solutions of p i,j,k+1 . At last, the calculation enters into the next time step t = t k+2 , and solutions w i,j,k+1 and p i,j,k+1 are regarded as the solutions of the previous time step.

Experiment Parameters
For numerical calculation the main parameters are given as follows: Young's modulus = 30 GPa, Poisson's ratio ν = 0.2, the injection rate at the wellbore Q 0 = 10 −4 m 3 /s, fluid viscosity µ = 10 −3 Pa·s. The initial conditions of the fracture are as follows: the initial radius R 0 = 0.5 m, the initial pressure p 0 = 0.1 MPa. Figure 4 shows the fixed grid of fracture, grid parameters ∆r = 0.02 m and ∆θ = 20 • (n = 360 • /20 • = 18). The symmetrical yellow areas represent one layer, i.e., the pay zone with a fracture toughness in the horizontal direction, the other white areas represent the other layers, i.e., the barrier layers with another fracture toughness. propagates, and the fracture planar area with the fracture elements is updated. Therefore, the entire calculation area for the model is updated, and the calculation needs to return to the first step for new values of , , , until all stress intensity factors ( ) satisfy the condition (18). The final values of , , as well as the new fracture front satisfying the conditions (11) and (18)  are inversely substituted into Equation (12) for the solutions of , , . At last, the calculation enters into the next time step = , and solutions , , and , , are regarded as the solutions of the previous time step.

Experiment Parameters
For numerical calculation the main parameters are given as follows: Young's modulus = 30 , Poisson's ratio = 0.2, the injection rate at the wellbore = 10 ⁄ , fluid viscosity = 10 ⋅ . The initial conditions of the fracture are as follows: the initial radius = 0.5 , the initial pressure = 0.1 . Figure 4 shows the fixed grid of fracture, grid parameters Δ = 0.02 and Δ = 20 ∘ ( = 360 ∘ 20 ∘ ⁄ = 18). The symmetrical yellow areas represent one layer, i.e., the pay zone with a fracture toughness in the horizontal direction, the other white areas represent the other layers, i.e., the barrier layers with another fracture toughness.  In this paper we want to study the effects of the fracture toughness contrast of the two different layers and the relative width of the pay zone, i.e., the proportion of the yellow areas on fracture propagation in early time, so the fracture toughness contrast and the proportion of the yellow areas are selected as the governing variables of numerical experiments. The crossover experiment strategy is adopted for quantitative analysis of the effects of these two governing variables.
In Figure 4 the proportion of the yellow areas involving only one pair of radial grids is ξ = 1/9. According to this grid discretization in Figure 4 we choose 3 values of the proportion: ξ = 1/9, 3/9 and 5/9 corresponding to 1, 3 and 5 pairs of radial grids as shown in Figure 5. The values of fracture toughness in the two layers for numerical experiments are obtained from papers [31,32]: the fracture toughness of barrier layers is fixed at K b IC = 2 MPa·m 0.5 , the fracture toughness of the pay zone changes from the minimum value 1 MPa·m 0. 5 Table 1.
In this paper we want to study the effects of the fracture toughness contrast of the two different layers and the relative width of the pay zone, i.e., the proportion of the yellow areas on fracture propagation in early time, so the fracture toughness contrast and the proportion of the yellow areas are selected as the governing variables of numerical experiments. The crossover experiment strategy is adopted for quantitative analysis of the effects of these two governing variables.
In Figure 4 the proportion of the yellow areas involving only one pair of radial grids is = 1 9 ⁄ . According to this grid discretization in Figure 4 we choose 3 values of the proportion: = 1 9 ⁄ , 3 9 ⁄ and 5 9 ⁄ corresponding to 1, 3 and 5 pairs of radial grids as shown in Figure 5. The values of fracture toughness in the two layers for numerical experiments are obtained from papers [31,32]: the fracture toughness of barrier layers is fixed at = 2 ⋅ . , the fracture toughness of the pay zone changes from the min-  The inhomogeneity of rock fracture toughness leads to the uneven propagation of fracture, i.e., different front points of the fracture may grow forward with different velocities, moreover, at each moment there are different front points, which satisfy the propagation condition and grow forward. In principle, the smaller the time step ∆t, the more accurately the simulation captures this uneven feature of fracture propagation. However, an exces-sively small-time step leads to non-convergence of model calculation, and an overly large time step makes it difficult to capture the uneven feature of fracture propagation. Therefore, in our experiments, the time step is assumed as = 0.01 s, and the total simulation time of each experiment is set as T = 20 s for the early time evolution of fracture propagation, i.e., the time step number: N = 2000.

Experiment Results
According to the above parameters, 27 numerical experiments are simulated based on this model and algorithm. Figures 6-8 show the results of the distribution of fluid pressure in the fracture at the moment = 20 s, Figure 9 shows comparison of fluid pressures along radial zones θ = 0 and θ = π/2 for different numerical experiments, and       The results of the fracture opening in Figures 10-12 indicate that the distribution of fracture opening is similar to the crack profile. The distributions of fracture opening along all radial zones are the same under the uniform propagation, see Figures 10 e, 11e and 12e. When the fracture has a very irregular contour, for example, Figures 10a and 12i, where the radius of the radial zone = 0 is much larger (or much smaller) than the radius of the radial zone = 2 ⁄ , the distributions of fracture opening along these two radial zones are very different from each other, see Figure 13a,f, because for all radial zones boundary conditions of the fracture opening at the fracture front ( = 0) and at the wellbore ( ) are the same (see the start and the end of all curves of fracture opening in Figure 13), so the average fracture opening gradient of the longer radial zone with a larger radius is smaller than the average fracture opening gradient of the shorter radial zone with a smaller radius. It means that because the radial zone with a less fracture toughness propagates faster than radial zones with a larger fracture toughness, the fracture opening gradient of the shorter radial zone near the fracture front deeply increases, which leads to the larger and larger stress intensity factor according to Equation (18) (see the red curves in Figure 13b,d,f). When the stress intensity factor of the shorter radial zone reaches the limit value ( ), it propagates forward and the radius increases. After that its stress intensity factor rapidly decreases, as shown by the black curves in Figure 13a,c,e. The mechanism of uneven fracture propagation is that the zones with less fracture toughness drive the zones with larger fracture toughness to propagate because of the compatibility of fracture opening on each radial zone. This is the reason why the radial zone with smaller fracture toughness cannot always propagate without the growth of the radial zone with larger fracture toughness, except under some special conditions. This feature is clearly seen in Figure 13a. Although the radial zone = 0 with the smaller fracture toughness propagated very fast from the initial radius 0.5 to 3.9 , the radial zone = 2 ⁄ with a larger fracture toughness did not keep still, and also propagated from 0.5 to 1.3 , just the velocity was smaller. It is opposite for the radial zones = 0 and = 2 ⁄ in Figure   13b: the radial zone = 0 propagated slower, while the radial zone = 2 ⁄ propagated faster.

Validation of Model and Algorithm
Savitski and Detournay [38] studied the propagation of a penny-shaped hydraulic fracture in the solid medium with a uniform fracture toughness, namely the experiment 5, 14 or 23 in Table 1. The large-toughness asymptotic solution to the fracture radius was obtained [38]: where the length scaling parameter and the dimensionless fracture radius are formulated:  Figure 14. It shows that the two results generally coincide with each other very well, with only a small difference at the beginning and the end of the time interval. The reasons are that in this model, the fracture has an initial radius of = 0.5 and the initial pressure distribution is not enough to propagate the fracture, so at the beginning of the simulation, the fracture radius remains unchanged. Furthermore, it should be noted that the asymptotic solution in [38] is obtained with the assumption that the fracture is always in dynamic propagation equilibrium; however, in our model, the numerical solution is obtained with the assumption that the fracture propagates radially one  Figure 9a shows the specific distribution of weak pressure zone on the radial zone θ = 0 near the fracture front. Its feature is a sharp decrease in the pressure gradient near the fracture front. The reason is that the fracture front located at these weak pressure zones propagated at the previous time step, which led to stress release and the emergence of new spaces near these fracture tips. The new spaces were initially liquid-free with zero pressure, and large pressure gradients between the propagating zones and the surrounding non-propagating areas were generated at these tips, i.e., these weak pressure zones. As the hydraulic fluid flows into these new spaces, the pressure gradients near the fracture front become smaller and smaller, as shown in Figure 9b-f. These features can only be observed in the uneven crack propagation. For the case of uniform crack propagation, as shown in Figures 6e, 7e and 8e, near the fracture front there is no distinct weak pressure zone, because the whole fracture tips at different locations propagate at the same time, and it is impossible to form distinct weak pressure area. It should be noted that these weak pressure areas only affect the pressure distribution near the fracture front and lead to the circular flow at these locations. Away from the fracture front, for example, near the wellbore, only the radial pressure gradient appears due to fluid injection, there is almost no pressure gradient along the hoop direction, as shown in Figures 6-8. This feature can be more clearly observed in Figure 9: pressures of different radial zones (θ = 0 and θ = π/2) near the coordinate origin are almost the same.
The results of the fracture opening in Figures 10-12 indicate that the distribution of fracture opening is similar to the crack profile. The distributions of fracture opening along all radial zones are the same under the uniform propagation, see Figure 10e, Figures 11e and 12e. When the fracture has a very irregular contour, for example, Figures 10a and 12i, where the radius of the radial zone θ = 0 is much larger (or much smaller) than the radius of the radial zone θ = π/2, the distributions of fracture opening along these two radial zones are very different from each other, see Figure 13a,f, because for all radial zones boundary conditions of the fracture opening at the fracture front (w = 0) and at the wellbore (w 0 ) are the same (see the start and the end of all curves of fracture opening in Figure 13), so the average fracture opening gradient of the longer radial zone with a larger radius is smaller than the average fracture opening gradient of the shorter radial zone with a smaller radius. It means that because the radial zone with a less fracture toughness propagates faster than radial zones with a larger fracture toughness, the fracture opening gradient of the shorter radial zone near the fracture front deeply increases, which leads to the larger and larger stress intensity factor according to Equation (18) (see the red curves in Figure 13b,d,f). When the stress intensity factor of the shorter radial zone reaches the limit value (K IC ), it propagates forward and the radius increases. After that its stress intensity factor rapidly decreases, as shown by the black curves in Figure 13a,c,e. The mechanism of uneven fracture propagation is that the zones with less fracture toughness drive the zones with larger fracture toughness to propagate because of the compatibility of fracture opening on each radial zone. This is the reason why the radial zone with smaller fracture toughness cannot always propagate without the growth of the radial zone with larger fracture toughness, except under some special conditions. This feature is clearly seen in Figure 13a. Although the radial zone θ = 0 with the smaller fracture toughness propagated very fast from the initial radius 0.5 m to 3.9 m, the radial zone θ = π/2 with a larger fracture toughness did not keep still, and also propagated from 0.5 m to 1.3 m, just the velocity was smaller. It is opposite for the radial zones θ = 0 and θ = π/2 in Figure 13b: the radial zone θ = 0 propagated slower, while the radial zone θ = π/2 propagated faster.

Validation of Model and Algorithm
Savitski and Detournay [38] studied the propagation of a penny-shaped hydraulic fracture in the solid medium with a uniform fracture toughness, namely the experiment 5, 14 or 23 in Table 1. The large-toughness asymptotic solution to the fracture radius was obtained [38]: where the length scaling parameter L k and the dimensionless fracture radius γ k are formulated: where K = 4(2/π) 1/2 K IC , γ k0 = 0.8546, γ k1 = −0.7349 and the dimensionless viscosity . The comparison of the results of the large-toughness asymptotic solution to the fracture radius and numerical experiments 5, 14 or 23 with the uniform rock toughness in this paper are exhibited in Figure 14. It shows that the two results generally coincide with each other very well, with only a small difference at the beginning and the end of the time interval. The reasons are that in this model, the fracture has an initial radius of R 0 = 0.5 m and the initial pressure distribution is not enough to propagate the fracture, so at the beginning of the simulation, the fracture radius remains unchanged. Furthermore, it should be noted that the asymptotic solution in [38] is obtained with the assumption that the fracture is always in dynamic propagation equilibrium; however, in our model, the numerical solution is obtained with the assumption that the fracture propagates radially one element at a time, so the difference between the asymptotic solution and the numerical solution increases with time, as shown in Figure 14 at the end of the time interval. solution increases with time, as shown in Figure 14 at the end of the time interval.
The reason for not choosing the assumption of dynamic propagation equilibrium for our model is that this assumption would involve one additional unknown for each boundary grid element into the model, namely the propagating length at each time step. However, there is no corresponding linearly independent equation introduced into the model for solving.  The fracture profile is generally a circle in Figure 15 under = 1 9 ⁄ . When the fracture toughness contrast < 1, the fracture profile is like a circle with a pair of obverse wedges at the symmetric pay zone with the smaller fracture toughness (see the red curve in Figure 15). As the fracture toughness contrast disappears ( → 1), the pair of wedges also gradually disappears, see the curves of = 0.625, 0.75 and 0.875. The green dotted curve in Figure 15 represents the standard circular fracture profile under the uniform propagation with no fracture toughness contrast ( = 1), which is the same as in Figures The reason for not choosing the assumption of dynamic propagation equilibrium for our model is that this assumption would involve one additional unknown for each boundary grid element into the model, namely the propagating length at each time step. However, there is no corresponding linearly independent equation introduced into the model for solving. ness (here we define a pair of reverse wedges with their tips facing each other, conversely it is a pair of obverse wedges). When = 1.25, the reverse wedges are not obvious, because the length of reverse wedges, namely, the retracted length of radial zones ( = 0 and = ) with the larger toughness is relatively small. With the increase of fracture toughness contrast ( = 1.5, 1.75 and 2.0), the shape of reverse wedges becomes more obvious, i.e., the retracted length increases with the increase of fracture toughness contrast under > 1. and even to 5 9 ⁄ , the fracture profile consisting of a circle with a pair of obverse wedges transforms into a horizontal rectangle with its long side along -axis. As the fracture toughness contrast disappears ( → 1), the difference between the long and short sides of this horizontal rectangle also disappears, and the fracture profile approaches a circle (the green dotted curve). With the increase of fracture toughness contrast over unity, the fracture profile is transformed into a vertical rectangle with the long side along -axis. The results of fracture profiles under = 3 9 ⁄ , 5 9 ⁄ in Figures 16 and 17 show that when the fracture toughness of the pay zone is larger than that of the barrier layers, the hydraulic fracture easily and mainly propagates toward the barrier layers, and the fracture along the pay zone is relatively small, which deviates from the goal of hydraulic fracturing.  The fracture profile is generally a circle in Figure 15 under ξ = 1/9. When the fracture toughness contrast ζ < 1, the fracture profile is like a circle with a pair of obverse wedges at the symmetric pay zone with the smaller fracture toughness (see the red curve in Figure 15). As the fracture toughness contrast disappears ( ζ → 1 ), the pair of wedges also gradually disappears, see the curves of ζ = 0.625, 0.75 and 0.875. The green dotted curve in Figure 15 represents the standard circular fracture profile under the uniform propagation with no fracture toughness contrast (ζ = 1), which is the same as in Figures 16 and 17. When the fracture toughness contrast ζ > 1, the fracture profile is like a circle losing a pair of reverse wedges at the symmetric pay zone with the larger fracture toughness (here we define a pair of reverse wedges with their tips facing each other, conversely it is a pair of obverse wedges). When ζ = 1.25, the reverse wedges are not obvious, because the length of reverse wedges, namely, the retracted length of radial zones (θ = 0 and θ = π) with the larger toughness is relatively small. With the increase of fracture toughness contrast (ζ = 1.5, 1.75 and 2.0), the shape of reverse wedges becomes more obvious, i.e., the retracted length increases with the increase of fracture toughness contrast under ζ > 1.

y (m)
The fracture profiles under the proportion of pay zone ξ = 3/9 and 5/9 are very similar to each other, but significantly different from the profiles under ξ = 1/9, see Figures 16 and 17. When ζ < 1, as the proportion of pay zone increases from 1/9 to 3/9, and even to 5/9, the fracture profile consisting of a circle with a pair of obverse wedges transforms into a horizontal rectangle with its long side along x-axis. As the fracture toughness contrast disappears ( ζ → 1 ), the difference between the long and short sides of this horizontal rectangle also disappears, and the fracture profile approaches a circle (the green dotted curve). With the increase of fracture toughness contrast over unity, the fracture profile is transformed into a vertical rectangle with the long side along y-axis. The results of fracture profiles under ξ = 3/9, 5/9 in Figures 16 and 17 show that when the fracture toughness of the pay zone is larger than that of the barrier layers, the hydraulic fracture easily and mainly propagates toward the barrier layers, and the fracture along the pay zone is relatively small, which deviates from the goal of hydraulic fracturing. In this model, the radii of radial zones = 0 and = 2 ⁄ are the most representative lengths to describe the fracture geometry because of symmetry, so here we take the radii of radial zone = 0 and = 2 ⁄ as the "half-length" ( ) and the "half-height" ( ) of the fracture according to the traditional practice in most models of hydraulic fracturing. The ratio of and is called the fracture shape coefficient = ⁄ . Figures 19-21 show the evolution of , and with time for the experiments of Groups 1, 2 and 3 corresponding to = 1 9 ⁄ , 3 9 ⁄ and 5 9 ⁄ , respectively. It is known from Figure 19 ( = 1 9 ⁄ ) that when < 1, is always greater than , but with the increase of the fracture toughness contrast from 0.5 to 1, the difference between and decreases and approaches zero. When > 1, is in turn greater than and the difference between and appears again, even rapidly increases with fracture toughness contrast changing from 1 to 2. The fracture shape coefficient initially changes rapidly and then approaches one certain limit value for all fracture toughness contrasts. The difference for < 1 and > 1 is that when < 1, is an increasing function of time from unit to the certain limit value, while > 1, is a decreasing function of time from unity to the limit value. According to the data in Figure 19, it can be estimated that the limit values of the fracture shape coefficient under =  By comparison, it is found that the fracture propagation is more sensitive to the formation with less than unity fracture toughness contrast: when the fracture toughness of pay zone is less and more than that of the barrier layers the same proportion, the excess propagation length on account of less than unity fracture toughness contrast is larger than the shortened length owing to more than unity fracture toughness contrast, for example, K p IC = 1 and 3.0 MPa·m 0.5 with a proportion of 50% less and more than the fracture toughness of barrier zones, the excess propagation lengths of the radial zone θ = 0 compared with the uniform propagation are 2.0, 1.8, 1.4 m corresponding to ξ = 1/9, 3/9 and 5/9, while the corresponding shortened lengths are only 0.5, 0.7, 0.8 m.
It is clearly observed in Figure 18 that the fracture toughness contrast between the pay zone and the barrier zones is the core factor of determining the overall propagation direction of the fracture. No matter what proportion the pay zone occupies (ξ = 1/9, 3/9 or 5/9), when the fracture toughness of the pay zone is less than that of the barrier zones (ζ < 1), the fracture in total propagates along the pay zone, see Figure 18a-d, on the contrary, when ζ > 1, the fracture propagates generally perpendicular to the pay zone, see Figure 18f-i. The proportion of pay zone only plays a role in strengthening or weakening the overall direction of fracture propagation, for example, in subgraphs (a-c) of Figure 18 the horizontal propagation of the fracture is very significant under ξ = 1/9, but only concentrated in very limited zones, the overall fracture profile is still a circle; when the Energies 2022, 15, 1967 21 of 35 proportion of pay zone is large enough (ξ = 3/9, 5/9), the fracture profiles are completely transformed from a circle into a horizontal rectangular. Furthermore, when the variables ζ and ξ are large enough, some fronts of the fracture cannot propagate, such as the blue curve of ζ = 2.0 and ξ = 5/9. It indicates that the combined effect of fracture toughness contrast and the proportion of the pay zone is very significant for the fracture profile.
It should be noted that the fracture profiles under ξ = 3/9, 5/9 in Figure 18a-c are very similar to the fracture shape of the classical Perkins-Kern-Nordgren (PKN) model [39,40]. This further proves the importance of the effects of the fracture toughness contrast and the proportion of pay zone on the fracture profile.
In this model, the radii of radial zones θ = 0 and θ = π/2 are the most representative lengths to describe the fracture geometry because of symmetry, so here we take the radii of radial zone θ = 0 and θ = π/2 as the "half-length" (R l ) and the "half-height" (R h ) of the fracture according to the traditional practice in most models of hydraulic fracturing. The ratio of R l and R h is called the fracture shape coefficient η = R l /R h . Figures 19-21 show the evolution of R l , R h and η with time t for the experiments of Groups 1, 2 and 3 corresponding to ξ = 1/9, 3/9 and 5/9, respectively. It is known from Figure 19 (ξ = 1/9) that when ζ < 1, R l is always greater than R h , but with the increase of the fracture toughness contrast from 0.5 to 1, the difference between R l and R h decreases and approaches zero. When ζ > 1, R h is in turn greater than R l and the difference between R l and R h appears again, even rapidly increases with fracture toughness contrast changing from 1 to 2. The fracture shape coefficient η initially changes rapidly and then approaches one certain limit value for all fracture toughness contrasts. The difference for ζ < 1 and ζ > 1 is that when ζ < 1, η is an increasing function of time from unit to the certain limit value, while ζ > 1, η is a decreasing function of time from unity to the limit value. According to the data in Figure 19, it can be estimated that the limit values of the fracture shape coefficient η under ξ = 1/9 are about 2.55, 1.58, 1.  Figure 19 under ξ = 1/9. The main difference is that with the increase of the proportion of pay zone, the velocities of fracture "half-length" and "half-height" change a lot, which greatly affects the evolution of the fracture shape coefficient. Table 2 gives the estimated limit value of the fracture shape coefficient with respect to the fracture toughness contrast and the proportion of the pay zone. From this table, we know that when ζ > 1, the limit values of η are inversely proportional to parameters ξ and ζ; when ζ < 1, the values of η are proportional only to the parameter ζ. Under the same fracture toughness contrast, the limit values of η reach the maximum with ξ = 3/9, not with ξ = 5/9. Comparing fracture profiles in Figures 15-18 with the limit values of fracture shape coefficient, we find that the formation with the moderate value of the proportion of the pay zone (ξ = 3/9) and the small fracture toughness contrast (ζ= 0.5) is the best for hydraulic fracturing, and the fracture profile is the most ideal.  ⁄ . The main difference is that with the increase of the proportion of pay zone, the velocities of fracture "half-length" and "half-height" change a lot, which greatly affects the evolution of the fracture shape coefficient. Table 2 gives the estimated limit value of the fracture shape coefficient with respect to the fracture toughness contrast and the proportion of the pay zone. From this table, we know that when > 1, the limit values of are inversely proportional to parameters and ; when < 1, the values of are proportional only to the parameter . Under the same fracture toughness contrast, the limit values of reach the maximum with = 3 9 ⁄ , not with = 5 9 ⁄ . Comparing fracture profiles in Figures 15-18 with the limit values of fracture shape coefficient, we find that the formation with the moderate value of the proportion of the pay zone ( = 3 9 ⁄ ) and the small fracture toughness contrast ( =0.5) is the best for hydraulic fracturing, and the fracture profile is the most ideal. The shape coefficient can be formulated as follows Figure 19. Relationships of the fracture "half-length" (R l ), "half-height" (R h ) and the fracture shape coefficient η with respect to time for Group 1 experiments with ξ = 1/9.
The shape coefficient η can be formulated as follows where v l and v h are the average propagation velocities of horizontal and vertical directions at the moment t. Another important feature of the curves of fracture shape coefficient η in Figures 19-21 is that: whether ζ < 1 or ζ > 1, the shape coefficient η always increases or decreases rapidly only at the initial period of fracture propagation, after which it changes (increases or decreases) very slightly and just fluctuates in a small range around the limit values. It means that the ratio between the horizontal and the vertical velocities changes rapidly at the beginning, and after approaching the limit value of η in Table 2, the ratio of velocities is basically unchanged. In a word, after the initial uneven propagation with a rapidly changed velocity rate, the fracture propagation will enter into a stable stage with an approximately constant velocity ratio between the horizontal and vertical directions, and in this stage, the fracture profile has the similarity with respect to time. Figures 22-24 give the evolution of fluid pressure (subgraphs (a)) and the fracture opening (subgraphs (b)) at the wellbore corresponding to experiments of Group 1, 2 and 3, in which ξ = 1/9, 3/9 and 5/9, respectively. In each figure of Figures 22-24 subgraphs (a) and (b) have the same legend, so it was drawn only in subgraph (a). around the limit values. It means that the ratio between the horizontal and the vertical velocities changes rapidly at the beginning, and after approaching the limit value of in Table 2, the ratio of velocities is basically unchanged. In a word, after the initial uneven propagation with a rapidly changed velocity rate, the fracture propagation will enter into a stable stage with an approximately constant velocity ratio between the horizontal and vertical directions, and in this stage, the fracture profile has the similarity with respect to time.  As we can see from Figures 22a, 23a and 24a, after entering into propagation stage pressure at the wellbore decreases with time because of fracture propagation. Furthermore, for all values ξ = 1/9, 3/9 and 5/9, the pressure at the wellbore p 0 increases with the increase of fracture toughness contrast. The fracture opening at the wellbore w 0 has the same evolution characteristic as the pressure at the wellbore, but there are two exceptions: ζ = 1.75 and ζ = 2.0, in which w 0 are smaller than that of experiments ζ = 1.5. These exceptions are evident in Figures 23 and 24, but not in Figure 22. The reason is that with the increase of the proportion of the pay zone (ξ = 3/9 and 5/9), the ability to hinder fracture propagation is significantly enhanced when ζ = 1.75 and ζ = 2.0 (see Figures 16 and 17: the pay zone with the larger fracture toughness propagated very slowly), even can lead to no propagation of the pay zone (see the fracture profile of experiment 27 with ξ = 5/9 and ζ = 2.0). Therefore, when the fracture toughness of the pay zone is large enough as ζ = 1.75 and ζ = 2.0, the difference between radii of the pay zone and the barrier zones is great, which leads to non-monotonically decrease of fracture opening with respect to the radial coordinate along some radial zones, even near the wellbore the fracture opening is a monotonically increasing function of the radial coordinate (see Figure 13c-f). This indicates that for some special cases (ξ = 3/9 and 5/9, ζ = 1.75 and 2.0) the fracture opening at the wellbore is not the maximum, which is significantly different from the 2D and 3D models for uniform fracture propagation.

Hydraulic Fracture Cleaning Process
The article considers the process of cleaning a fracture from hydraulic fracturing fluid. This process is performed by oil under the effect of inverse pressure drop, and is regarded within the frame of a planar 2-D approximation in the plane orthogonal to a vertical fracture [41]. The following geometry of the position of the wells is considered: four injection wells in the corners of the computational domain and one production well in the center. The central well is connected with a symmetrical fracture filled in with proppant. The reservoir area around the crack is filled with fracturing fluid. Using numerical modeling on the basis of the constructed mathematical model, a relationship is established between the quality of hydraulic fracture cleaning and the geometrical parameters of the fracture and the region filled with the hydraulic fracturing fluid.

Mathematical Model for Seepage Flows of Multi-Phase Fluids
The problem was solved under the following assumptions: 1.
The viscosity of the displaced fluid is greater than the viscosity of the displacing fluid in the reservoir; 3.
Capillary effects at the interface are not taken into account; 4.
Permeability and porosity are location dependent; 5.
Thermal effects and gravity are not taken into account; 6.
The outer boundary is impermeable; 7.
Seepage is modeled by Darcy's law, taking into account the relative permeabilities of the phases 8.
Hydraulic fracture-area of increased permeability and porosity 9.
Permeability has a random "ripple" that contributes to the onset of displacement instability.
The following system of equations was solved. For each fluid phase, the law of mass conservation is written as: Here, ϕ is porosity, s is saturation, is density (intrinsic), u k,j is j-th component of the seepage velocity of the k-th phase.
The average volumetric seepage velocity of the fluid is defined as: Summation of Equations (21), taking into account the liquid's incompressibility and constant in time porosity, leads to: Darcy's law for each phase is as follows: where µ k is the dynamic viscosity of the phase, K is absolute permeability of the medium, K R k is the relative permeability of k-th phase, p is the pressure in the pores. Substitution of Equations (24) into (22) leads to: Energies 2022, 15, 1967 27 of 35 where m k = k r k µ k -relative mobility. Substituting expression (24) into equation (23) gives an equation for the pressure: Equation (21), taking into account the introduced definitions (22,25,26), leads to an equation for saturation dynamics: where f k = m k m 1 +m 2 is the proportion of relative phase mobility. The relative permeabilities K R k are calculated using the Brooks-Corey model [42]: Here, k 0 k > 0 and n 0 k > 0 are the model parameters, and the effective saturation S k is determined by the residual saturations of 0 ≤ s res k ≤ 1 (s res 1 + s res 2 < 1). A hydraulic fracture is modeled as an area of increased permeability and porosity. Thus, absolute permeability and porosity are expressed as: where F is the set of points lying inside the fracture. It is also worth noting that a weak random "ripple" is superimposed on the permeability, so that finally the absolute permeability is defined as: where ξ is a random variable uniformly distributed in the interval [−1; +1], δ is a rather small quantity. This random "ripple" contributes to the onset of displacement instability. Figure 25 shows the geometry of the model problem. The reservoir is initially saturated with pore fluid (oil). There are four wells in the corners with constant pressure in it. In the center of the area there is a production well and a hydraulic fracture filled with proppant and hydraulic fracturing fluid. It is assumed that part of the fracturing fluid has leaked into the surrounding rock. It is also assumed here that the region impregnated with the hydraulic fracturing fluid is elliptical (region A). The problem under consideration is cleaning the region A from the remains of hydraulic fracturing fluid.
Due to the pressure drop, pore fluid is filtered in an isotropic porous medium to the producing well and to the outside. The viscosity of the fracturing fluid is greater than the viscosity of the oil. Therefore, the process is unstable. The outer walls are impermeable.
The boundary conditions: Σ in : s = 0, P = P in Σ out : P = P out < P in Σ w : u n = 0 → ∂p ∂n Г = 0 (32) Energies 2022, 15, 1967 28 of 35 The initial conditions: The model consists of Equations (25)- (28), with the conditions (29)- (33).  Figure 25 shows the geometry of the model problem. The reservoir is initially saturated with pore fluid (oil). There are four wells in the corners with constant pressure in it. In the center of the area there is a production well and a hydraulic fracture filled with proppant and hydraulic fracturing fluid. It is assumed that part of the fracturing fluid has leaked into the surrounding rock. It is also assumed here that the region impregnated with the hydraulic fracturing fluid is elliptical (region A). The problem under consideration is cleaning the region A from the remains of hydraulic fracturing fluid. Due to the pressure drop, pore fluid is filtered in an isotropic porous medium to the producing well and to the outside. The viscosity of the fracturing fluid is greater than the viscosity of the oil. Therefore, the process is unstable. The outer walls are impermeable.

Effect of Fracture Length on Displacement Dynamics
The purpose of the first numerical experiment is to investigate the effect of fracture length ( ) on the cleaning process dynamics. Two calculations were carried out with different fracture lengths. The calculation parameters are given in Table 3.

Effect of Fracture Length on Displacement Dynamics
The purpose of the first numerical experiment is to investigate the effect of fracture length (L f racture ) on the cleaning process dynamics. Two calculations were carried out with different fracture lengths. The calculation parameters are given in Table 3. The dependence in Figure 26 initially behaves linearly, after that it continues to decrease, but never reaches a constant: the asymptote in this case is s = 0. Therefore, any sufficiently small or satisfactory saturation value can be taken as a criterion when it is considered that the fracture is already cleaned. If we assume that the fracture is already cleaned when the saturation of the hydraulic fracturing fluid reached 0.2, then it took about 5 days 8 h to clean the fracture 20 m long, and it took 9 days 13 h to clean the fracture 30 m long. It turns out that the ratio is 1.79. Therefore, the length of the fracture significantly affects the duration of the cleaning process. Figures 27-29 show the distribution patterns of the oil seepage velocity for various moments of time. The instability of displacement is observed in all figures.
In Figure 27, the largest flows are observed in the fracture, as well as near the production well. It can be seen that for a larger fracture, the fluid flows faster in the area near the well.
In Figure 28, for the case when the length is 20 m, the flow occurs along the entire fracture. In the case of a longer fracture, a breakthrough occurred in the area of the production well. Thus, at this time moment, not the entire hydraulic fracture is involved in the seepage process. In Figure 29 it can be noted that seepage occurs in almost the entire area for both fracture lengths. Oil is accumulated strongly in the rock near the well, continues to flow into it. Besides, there are areas that have not changed much: the accumulation of oil and the displacement of hydraulic fracturing fluid in them occurs very slowly.    The dependence in Figure 26 initially behaves linearly, after that it continues to decrease, but never reaches a constant: the asymptote in this case is s = 0. Therefore, any sufficiently small or satisfactory saturation value can be taken as a criterion when it is considered that the fracture is already cleaned. If we assume that the fracture is already cleaned when the saturation of the hydraulic fracturing fluid reached 0.2, then it took about 5 days 8 h to clean the fracture 20 m long, and it took 9 days 13 h to clean the fracture  In Figure 27, the largest flows are observed in the fracture, as well as near the production well. It can be seen that for a larger fracture, the fluid flows faster in the area near the well.     In Figure 27, the largest flows are observed in the fracture, as well as near the production well. It can be seen that for a larger fracture, the fluid flows faster in the area near the well.    In Figure 27, the largest flows are observed in the fracture, as well as near the production well. It can be seen that for a larger fracture, the fluid flows faster in the area near the well.

Effect of Fracture Fluid Viscosity on Displacement Dynamics
The purpose of the second numerical experiment is to investigate the effect of fracture fluid viscosity (µ2) on the cleaning process dynamics. Two calculations were carried out with different viscosities of fracturing fluid. The calculation parameters are given in Table 4. In Figure 30 we can see how the saturation of the displaced fluid in the elliptical region changes with time. It is noticeable that on each graph three parts can be distinguished. The first part, similar to the linear part, ends when there is a breakthrough of oil to the production well. The second part is the curved up section, and the third is the curved down section.
In Figure 28, for the case when the length is 20 m, the flow occurs along the entire fracture. In the case of a longer fracture, a breakthrough occurred in the area of the production well. Thus, at this time moment, not the entire hydraulic fracture is involved in the seepage process. In Figure 29 it can be noted that seepage occurs in almost the entire area for both fracture lengths. Oil is accumulated strongly in the rock near the well, continues to flow into it. Besides, there are areas that have not changed much: the accumulation of oil and the displacement of hydraulic fracturing fluid in them occurs very slowly.

Effect of Fracture Fluid Viscosity on Displacement Dynamics
The purpose of the second numerical experiment is to investigate the effect of fracture fluid viscosity (µ2) on the cleaning process dynamics. Two calculations were carried out with different viscosities of fracturing fluid. The calculation parameters are given in Table 4. In Figure 30 we can see how the saturation of the displaced fluid in the elliptical region changes with time. It is noticeable that on each graph three parts can be distinguished. The first part, similar to the linear part, ends when there is a breakthrough of oil to the production well. The second part is the curved up section, and the third is the curved down section.  The section curved up is the time interval when the displacement velocity increases. This phenomenon can be explained by the fact that during this time interval, the liquid is displaced from those areas that were not previously involved. These areas are shown in Figure 31. Then, when these hard-to-reach areas are partially cleared, displacement again slows down (third area). For the considered cases, the difference between the two graphs ( Figure 30) is not as significant as for the case of different fracture lengths. Using the same criterion (s = 0.2), it took about 5 days 8 h to clean the fracture for µ 2 = 0.1 Pa s, and it took 4 days 11 h to clean the fracture for µ 2 = 0.08 Pa s. It turns out that the ratio is 0.84, which is greater than a viscosity ratio of 0.8. Therefore, fracturing fluid viscosity has a slightly lesser effect on the duration of the cleaning process, than fracture length.
The section curved up is the time interval when the displacement velocity increases. This phenomenon can be explained by the fact that during this time interval, the liquid is displaced from those areas that were not previously involved. These areas are shown in Figure 31. Then, when these hard-to-reach areas are partially cleared, displacement again slows down (third area). For the considered cases, the difference between the two graphs ( Figure 30) is not as significant as for the case of different fracture lengths. Using the same criterion (s = 0.2), it took about 5 days 8 h to clean the fracture for = 0.1 • s, and it took 4 days 11 h to clean the fracture for = 0.08 • s. It turns out that the ratio is 0.84, which is greater than a viscosity ratio of 0.8. Therefore, fracturing fluid viscosity has a slightly lesser effect on the duration of the cleaning process, than fracture length.  Figure 32 shows the presence of a breakthrough in the case of the high viscosity and the absence of a breakthrough in the case of the low viscosity. It turns out that a breakthrough for higher viscosity occurs earlier than for lower viscosity. When the fracturing fluid has a high viscosity, a breakthrough through areas near the production well is possible. Thus, the hydraulic fracture effectively does not work, since the fracture is not connected to the seepage process.   Figure 32 shows the presence of a breakthrough in the case of the high viscosity and the absence of a breakthrough in the case of the low viscosity. It turns out that a breakthrough for higher viscosity occurs earlier than for lower viscosity. When the fracturing fluid has a high viscosity, a breakthrough through areas near the production well is possible. Thus, the hydraulic fracture effectively does not work, since the fracture is not connected to the seepage process.
The section curved up is the time interval when the displacement velocity increases. This phenomenon can be explained by the fact that during this time interval, the liquid is displaced from those areas that were not previously involved. These areas are shown in Figure 31. Then, when these hard-to-reach areas are partially cleared, displacement again slows down (third area). For the considered cases, the difference between the two graphs ( Figure 30) is not as significant as for the case of different fracture lengths. Using the same criterion (s = 0.2), it took about 5 days 8 h to clean the fracture for = 0.1 • s, and it took 4 days 11 h to clean the fracture for = 0.08 • s. It turns out that the ratio is 0.84, which is greater than a viscosity ratio of 0.8. Therefore, fracturing fluid viscosity has a slightly lesser effect on the duration of the cleaning process, than fracture length.  Figure 32 shows the presence of a breakthrough in the case of the high viscosity and the absence of a breakthrough in the case of the low viscosity. It turns out that a breakthrough for higher viscosity occurs earlier than for lower viscosity. When the fracturing fluid has a high viscosity, a breakthrough through areas near the production well is possible. Thus, the hydraulic fracture effectively does not work, since the fracture is not connected to the seepage process.

Discussion
A planar-3D non-uniform propagation model with consideration of rock's inhomogeneous fracture toughness has been developed for studying the effect of inhomogeneous fracture toughness on the fracture propagation. This model allows fluid flow along the radial and angular directions and considers the uneven propagation at different fracture tip locations. The inhomogeneous fracture toughness of formation is described by two parameters (the fracture toughness contrast ζ and the proportion of the pay zone ξ). According to these two controlling parameters, three groups of numerical experiments are designed and classified by the proportion of pay zone, in each group, there are nine subgroups of experiments corresponding to nine different levels of fracture toughness contrast. Computer codes for these numerical experiments were developed and implemented using the Matlab software package. The code was verified by comparing with exact solutions of a 2-D plane fracture problem.
Numerical results of experiments show that the fracture propagation properties including the fracture geometry, fluid pressure distribution and the fracture opening are significantly dependent on these two controlling parameters, especially on the fracture toughness contrast, which mainly controls the overall direction of fracture propagation: when ζ < 1, i.e., the fracture toughness of pay zone is less than that of barrier zones, the fracture mainly propagates along the pay zone (along the horizontal direction, R l > R h ); when ζ > 1, i.e., the fracture toughness of pay zone is larger than that of barrier zones, the fracture mainly propagates along the barrier zones (along the vertical direction, R l < R h ). Besides, the controlling parameter ξ mainly enhances the capability of the fracture toughness contrast determining the overall fracture direction: when the value of ξ is small (ξ = 1/9), no matter what the fracture toughness contrast is, it can only affect the partial fracture profile, the overall fracture profile is still a circle with or without a pair of narrow wedges. When ξ increases to 3/9, or 5/9, the combined effect of parameters ξ and ζ on the fracture profile is very significant: the fracture profile is transformed into an approximate ellipse or an approximate rectangle under some values of ζ (see Figures 16 and 17). The comparison of the results of the fracture "half-length" (R l ) and the fracture "half-height" (R h ) indicates that the velocity difference of uneven fracture propagation between the horizontal and vertical directions is mainly accumulated in the early stage of fracture propagation, after which the fracture non-uniformly propagates with an approximately constant velocity ratio between the horizontal and vertical directions. The analysis of the fluid pressure and the fracture opening at the wellbore shows that the greater the fracture toughness contrast, the larger the fluid pressure at the wellbore is obtained; the fracture opening has the similar evolution characteristic with two exceptions ζ = 1.75 and ζ = 2.0, because under the two values of fracture toughness contrast, the pay zone almost does not propagate, which leads to partially monotonous increase of fracture opening near the wellbore along the radius (see Figure 13), which further leads to some decrease of fracture opening at the wellbore. It should be noted that the presented results show the dependence of the hydraulic fracturing process only on the inhomogeneity of the fracture toughness of the medium, while in real conditions such inhomogeneity will also be associated with other inhomogeneities, such as stiffness or permeability and porosity. On the other hand, the obtained results show that the inhomogeneity of fracture toughness should not be ignored when taking into account real heterogeneous environment without proper justification.
In the hydraulic fracture cleaning process, when the fracturing fluid is replaced by oil from the deposit, instability of the displacement front takes place in the case of liquid with higher viscosity (fracturing fluid) being displaced by a liquid of lower viscosity (oil). In a small fracture, the cleaning process is faster, and the seepage process looks uniform. When the fracturing fluid has a higher viscosity than the displacing agent does, it is possible that the breakthrough will happen near the well, and not through the fracture. That means that the hydraulic fracture does not work properly as an oil collector; the fracture remains filled in by fracturing fluid not being involved in the oil seepage process. The breakthrough of the fluid to the production well is the faster the greater is the viscosity of the fracturing fluid. For the investigated hydraulic fracturing fluids with viscosities µ = 100 Mpa * s and µ = 80 Mpa * s, the effect of fracture length on the average saturation of hydraulic fracturing fluid is more noticeable than the effect of fracturing fluid viscosity. Due to the simplified problem statement, quantitative results should be treated with caution, but qualitative results should be kept in mind when evaluating the effectiveness of hydraulic fracturing.