Numerical Study of Elasto-Plastic Hydraulic Fracture Propagation in Deep Reservoirs Using a Hybrid EDFM–XFEM Method

: Rock yielding may well take place during hydraulic fracturing in deep reservoirs. The prevailing models based on the linear elastic fracture mechanics (LEFM) are incapable of describing the evolution process of hydraulic fractures accurately. In this paper, a hydro-elasto-plastic model is proposed to investigate the hydraulic fracture propagation in deep reservoirs. The Drucker–Prager plasticity model, Darcy’s law, cubic law and cohesive zone model are employed to describe the plastic deformation, matrix ﬂow, fracture ﬂow and evolution of hydraulic fractures, respectively. Combining the embedded discrete fracture model (EDFM), extended ﬁnite element method (XFEM) and ﬁnite volume method, a hybrid numerical scheme is presented to carry out simulations. A dual-layer iterative procedure is developed based on the ﬁxed-stress split method, Picard iterative method and Newton–Raphson iterative method. The iterative procedure is used to deal with the coupling between nonlinear deformation with fracture extension and ﬂuid ﬂow. The proposed model is veriﬁed against analytical solutions and other numerical simulation results. A series of numerical cases are performed to investigate the inﬂuences of rock plasticity, internal friction angle, dilatancy angle and permeability on hydraulic fracture propagation. Finally, the proposed model is extended to simulate multiple hydraulic fracture propagation. The result shows that plastic deformation can enhance the stress-shadowing effect. model of this problem includes the following three components: (i) rock deformation, (ii) ﬂuid ﬂow, and (iii) propagation criterion.


Introduction
With the rapid growth in the global demand for oil and gas resources and the successful development of middle-shallow formations (<4500 m), the exploitation of oil and gas reservoirs is gradually extended to the deep formations (>4500 m) [1,2]. Hydraulic fracturing technology, which can effectively improve the permeability of reservoir rocks, plays an important role in the development of oil and gas reservoirs [3]. It is important to predict the shape and trend of hydraulic fractures accurately. However, as deep reservoirs are characterized by high pressure, high temperature and high stress, the rock brittleness decreases and nonlinear deformation becomes obvious [4,5]. The prevailing fracture propagation models based on the linear elastic fracture mechanics are no longer applicable for predicting the evolution process of hydraulic fractures in deep ductile formations. Accounting for rock plasticity, establishing an efficient model to perform accurate simulation of hydraulic fracture propagation is significant for the development of deep oil and gas resources.
In the past decades, quite a few research works on hydraulic fracture propagation have been conducted, which can be roughly divided into three stages. In the first stage (1960s-1970s), researchers mainly focused on analytical models and solutions for the twodimensional hydraulic fracture, such as the classical KGD model [6,7] and PKN [8,9] model. In the second stage (1980s-1990s), on the basis of 2D fracture models, many quasi-3D and true-3D models [10,11] were presented. In the third stage (since the 1990s), with increasing complexity of problems and the development of computer technology, numerical models gradually replaced analytical models for the study of fluid-driven fracture propagation due to the limitations of analytical models. For hydraulic fracturing simulation, a key issue is to solve the discontinuous fields, including the displacement and pressure fields. Numerical methods usually used to simulate the fields with discontinuities are as follows: the finite element method [12,13], extended finite element method [14], discontinuous displacement method [15,16], distinct element method [17,18], discrete element model [19], lattice element model [20], discrete fracture network [21,22], and embedded discrete fracture model [23]. Their advantages and disadvantages are detailed in the literatures [24,25]. In the above methods, the XFEM allows fracture growth across the grid without remeshing, so it is widely used to solve discontinuous displacement fields. In addition, the EDFM, in which the fracture is embedded in non-matching grids, is usually employed to model the fluid flow between the fractures and the matrix in fractured formation. In this paper, the XFEM and EDFM are adopted to model the discontinuous displacement and pressure fields, respectively. Combined with the common advantages of the two, the limitations of the mesh size and shape on hydraulic fracturing simulation are relieved.
In most available fracture models, the rock is assumed to be linear elastic media. Few research works focus on hydraulic fracture propagation in elasto-plastic media. On the basis of the Mohr-Coulomb yielding theory, Papanastasiou [26] built an efficient cohesive finite element model to simulate hydraulic fracture propagation in inelastic media. Combining the extended finite element method and the cohesive zone model, Liu [27] developed a stabilized model to simulate hydraulic fracture propagation in brittle and ductile formations. Zeng [28] studied the effects of plastic deformation on fluid-driven fracture extension with the Drucker-Prager yield criterion. Liu [23] presented an efficient numerical model for elasto-plastic hydraulic fracture extension considering thermal effects. However, in these studies, fluid flow in the matrix is ignored by assuming an impermeable rock material, which has been proven to be an unrealistic assumption by considerable field evidence [29]. Sarris [30] added the seepage flow to Papanastasiou's model. Wang [31] carried out poro-elastic and poro-plastic modeling of hydraulic fracturing in brittle and ductile formations with the finite element method (FEM). In both studies of Sarris and Wang, a fully coupling solution method was adopted to conduct numerical simulations, which needed huge computational costs. Besides, in those studies mentioned above, only the simulation of straight hydraulic fracture extension was studied. In this work, the proposed model can simulate planar or non-planar elasto-plastic hydraulic fracture propagation in non-matching grids without remeshing.
In our previous studies, combining the XFEM and EDFM, multiphase flow in deformational fractured shale reservoirs was simulated [32,33]. Then, the hybrid EDFM-XFEM method was proposed to build an efficient hydro-mechanical model of hydraulic fracture propagation in brittle [34] and anisotropic [35] formations. In this work, to study fracture propagation in deep reservoirs, we extend the model to ductile formations accounting for the rock plasticity, and construct a dual-layer iterative procedure for solving the nonlinear problem caused by inelasticity and fracture propagation.

Mathematical Model
We consider the problem that a hydraulic fracture propagates in the deep reservoir, as shown in Figure 1. The reservoir is assumed as an isotropic and homogeneous porous media. Moreover, plastic deformation of the rock is taken into account. The mathematical We consider the problem that a hydraulic fracture propagates in the deep r as shown in Figure 1. The reservoir is assumed as an isotropic and homogeneou media. Moreover, plastic deformation of the rock is taken into account. The math model of this problem includes the following three components: (i) rock deform fluid flow, and (iii) propagation criterion.

Rock Deformation
The governing equations of rock deformation include the equilibrium equat metric equation and constitutive equation. The equilibrium equation is written as where σ is the total stress tensor and f is the body force vector.
According to Biot's theory, effective stress can be expressed as the following where σ' is the effective stress; pm is the fluid pressure in the matrix; and α is coefficient.
The geometric equation is written as follows: ( ) where ε, and u denote the strain and displacement, respectively. For the description of plastic deformation, the constitutive equation is chara by the Drucker-Prager plastic theory. The yielding function is written as follows where I1 and J2 represent the first invariant of the stress tensor and the second inv the deviatoric stress tensor, respectively, and a and k are the constants, express following: in which c and φ are the rock cohesion and internal friction angle, respectively. Based on the non-associated flow rule, the constitutive relationship can be e in the incremental form as follows:

Rock Deformation
The governing equations of rock deformation include the equilibrium equation, geometric equation and constitutive equation. The equilibrium equation is written as follows: where σ is the total stress tensor and f is the body force vector. According to Biot's theory, effective stress can be expressed as the following [36]: where σ is the effective stress; p m is the fluid pressure in the matrix; and α is the Biot coefficient. The geometric equation is written as follows: where ε, and u denote the strain and displacement, respectively. For the description of plastic deformation, the constitutive equation is characterized by the Drucker-Prager plastic theory. The yielding function is written as follows [37]: where I 1 and J 2 represent the first invariant of the stress tensor and the second invariant of the deviatoric stress tensor, respectively, and a and k are the constants, expressed as the following: in which c and ϕ are the rock cohesion and internal friction angle, respectively. Based on the non-associated flow rule, the constitutive relationship can be expressed in the incremental form as follows: where D and D ep represent the elastic and elasto-plastic stiffness tensor, respectively; Q denotes the flow potential; and A denotes the hardening parameter, which is set to be 0 in this paper. More details about the Drucker-Prager yielding criterion and the derivation of elasto-plastic constitutive relation can be found in our previous paper [23].

Fluid Flow
Fluid flow can be divided into the following three parts: flow in the matrix, flow in the fracture and the leak-off between them. The EDFM is adopted to characterize the fluid flow process. For matrix flow, the mass conservation equation can be written as follows: where ε v is the volumetric strain; M is Biot's modulus; q mf is the leak-off flow; S m is the area of the matrix gird; and v m is the fluid flow velocity described by Darcy's law, expressed by the following: in which k m and µ denote the matrix permeability and fluid viscosity, respectively. For fluid flow in the fracture, the mass conservation equation with consideration of the fracture width variation can be written as follows: where l mf is the length of the fracture grid; w f is the width of the fracture; and v f denotes the average fracture flow velocity, which can be given as the following: In addition, based on the EDFM, the leak-off flow can be calculated by the following: where d mf denotes the equivalent distance from the matrix element to the fracture segment [32], expressed as the following: in which x n represents the vertical distance from any one point in the matrix element to the fracture segment.

Cohesive Zone Model
Due to the plastic deformation of deep rocks, the fracture propagation criteria based on the LEFM, including the maximum circumferential stress criteria and J integrals, are no longer applicable for deep formations. The cohesive zone model is considered to be the most reliable propagation criterion in nonlinear mechanics [24,38], and it can be easily used with the FEM [39] and XFEM [40]. Here, the cohesive zone model is employed to describe the process of fracture propagation. In the cohesive zone model, the hydraulic fracture is decomposed into the following three zones: (i) non-damage zone, (ii) fracture process zone, and (iii) broken zone, as illustrated in Figure 2a. The bilinear traction-separation model is employed to describe the relationship between the traction and displacement jumps across the fracture. As shown in Figure 2b, t c is the cohesive traction in the following: t c ={t n , t s } T , in which t n and t s refer to the normal and shear components, respectively;u is the separation in the following: u = u + − u − . The area under the traction-separation curve Energies 2021, 14, 2610 5 of 18 equals the critical strain energy rate G C . The constitutive law of the cohesive zone model can be expressed as the following [41]: (13) where T is the tangential stiffness matrix of the traction-separation law.
process zone, and (iii) broken zone, as illustrated in Figure 2a. The bilinear traction-sep ration model is employed to describe the relationship between the traction and displac ment jumps across the fracture. As shown in Figure 2b, tc is the cohesive traction in th following: tc ={tn, ts} T , in which tn and ts refer to the normal and shear components, respe tively;   u is the separation in the following:   + − = − u u u . The area under the traction separation curve equals the critical strain energy rate GC. The constitutive law of the co hesive zone model can be expressed as the following [41]: where T is the tangential stiffness matrix of the traction-separation law.

Numerical Implementation
In this section, the mathematical model is discretized based on the XFEM, EDFM an FVM. Besides, combining the fixed-stress split method, Picard iterative method and New ton-Raphson iterative method, a dual-layer iterative procedure is established to solve th strongly nonlinear discretized equations. Based on the proposed numerical scheme, th numerical simulations in the subsequent part are conducted by the mathematics softwar MATLAB.

Displacement Field
The XFEM is adopted to simulate the inelastic stress field with the discontinuit based on the incremental theory. According to the XFEM, at each loading step, the expre sion of the incremental displacement field is written as follows: Δ Δ  denote the standard and enriched nodal displacement increments, r spectively; N is the shape function; and H(*) is the Heaviside function. For writing con venience, the Equation (14) can be simplified as follows: Accordingly, the incremental strain vector can be given as follows: where B and B enr denote the derivatives of the shape function matrices N and N enr , respe tively.
Using the Galerkin method, the weak form of Equation (1) can be obtained as follow

Numerical Implementation
In this section, the mathematical model is discretized based on the XFEM, EDFM and FVM. Besides, combining the fixed-stress split method, Picard iterative method and Newton-Raphson iterative method, a dual-layer iterative procedure is established to solve the strongly nonlinear discretized equations. Based on the proposed numerical scheme, the numerical simulations in the subsequent part are conducted by the mathematics software MATLAB.

Displacement Field
The XFEM is adopted to simulate the inelastic stress field with the discontinuity based on the incremental theory. According to the XFEM, at each loading step, the expression of the incremental displacement field is written as follows: where ∆u i , ∆ u j denote the standard and enriched nodal displacement increments, respectively; N is the shape function; and H(*) is the Heaviside function. For writing convenience, the Equation (14) can be simplified as follows: Accordingly, the incremental strain vector can be given as follows: where B and B enr denote the derivatives of the shape function matrices N and N enr , respectively. Using the Galerkin method, the weak form of Equation (1) can be obtained as follows: Based on the incremental method, the following applies: where subscript k denotes the loading step. Substituting Equation (18) into (17), the Newton-Raphson iterative form of the discretized system for solving incremental displacements can be obtained as follows: in which K and Q denote the stiffness matrix and coupling matrix, respectively; ∆F, F ext and F int denote, respectively, the incremental, external and internal force vector, expressed as the following:

Pressure Field
The finite volume method (FVM) is employed to simulate the discontinuous pressure field based on the EDFM. According to FVM, the mass conservation Equation (7) is integrated over the control volume and time, expressed as the following: Applying the Gaussian divergence theorem on the first term, and using the fully implicit scheme to conduct temporal discretization, Equation (21) becomes the following: where superscript n denotes the time step. Substituting Equations (8) and (11) into (22), the discrete equation can be given as the following:start a new page without indent 4.6cm in which the following applies: Similarly, the discrete form of fluid flow in the fracture can be obtained as follows: in which the following applies: Therefore, the final discrete scheme can be written in the form of the matrix as follows: in which the following applies:

Iterative Procedure
The nonlinearity of the established model is manifested in the following two aspects: (1) hydro-mechanical coupling and (2) the rock inelastic stress-strain relationship, which means that the solution involves two levels of the iterative process. Accordingly, a duallayer iterative procedure is built up for solving this strong nonlinear problem.
In the first layer, combining the Picard iterative method and fixed-stress split method, a global iterative algorithm is constructed to deal with nonlinear problems caused by the coupling of rock deformation and fluid flow, as shown in Figure 3. The fixed-stress split method was proved to be unconditionally stable [42] and robust for highly nonlinear problems [43]. Here, it is used to deal with a hydro-mechanical coupling problem. In accordance with the fixed-stress split method, volume strain is written as follows: Substituting Equation (29) into (25), after some manipulations, the fixed-stress iterative scheme of conservation in Equation (7) can be given as follows: Energies 2021, 14, 2610 8 of 18 in which the following applies: ; The Picard iterative method [29] is employed to handle the coupling problem between the fracture flow and varying width at every time step. If the convergence criterion is not satisfied, the fracture width is modified at each iteration step. The fluid pressure is recalculated with the modified fracture width from the previous iteration step. The convergence criterion is checked again after computing the new displacement based on the fluid pressure. Thus, the final discrete scheme (Equation (27)) can be derived as follows: in which the following applies: The Picard iterative method [29] is employed to handle the coupling problem between the fracture flow and varying width at every time step. If the convergence criterion is not satisfied, the fracture width is modified at each iteration step. The fluid pressure is recalculated with the modified fracture width from the previous iteration step. The convergence criterion is checked again after computing the new displacement based on the fluid pressure.
In the second layer, as illustrated in Figure 4, the nonlinear stress field Equation (19) is solved by the Newton-Raphson iterative algorithm, in which the stress state is updated by an implicit elastic predictor and a plastic corrector return-mapping scheme [37,44]. , x FOR PEER REVIEW 9 of 26 In the second layer, as illustrated in Figure 4, the nonlinear stress field Equation (19) is solved by the Newton-Raphson iterative algorithm, in which the stress state is updated by an implicit elastic predictor and a plastic corrector return-mapping scheme [37,44].

Numerical Results and Discussion
In this section, the proposed model is first validated through two cases of hydraulic fracture propagation. Then, the influences of rock plastic deformation and permeability on hydraulic fracture propagation in deep reservoirs are discussed. Finally, the proposed model is extended to simulate the simultaneous propagation of multiple hydraulic fractures.

Model Verification
Since no analytical solution of elasto-plastic hydraulic fracture propagation can be referred to, the correctness and effectiveness of the presented model and algorithm are verified from the following two aspects: (1) fracture propagation under the elastic deformation-KGD model, and (2) fracture propagation under plastic deformation-ABAQUS.

Fracture Propagation under Elastic Deformation
Neglecting the plastic deformation and matrix flow, the proposed model is simplified to the classical KGD problem [6,7]. The numerical results are now compared with the existing analytical solutions under elastic deformation. The input parameters are set as follows: E = 25 GPa, ν = 0.2, tc = 0.5 MPa, GC = 50 N/m, Qinject = 0.0005 m 3 /s, and μ = 0.1 Pa·s. Figure 5 shows the comparisons of the fracture half-length and the fracture width at wellbore between the numerical solutions and analytical solutions [45]. It can be found that the numerical results agree well with the analytical solutions.

Numerical Results and Discussion
In this section, the proposed model is first validated through two cases of hydraulic fracture propagation. Then, the influences of rock plastic deformation and permeability on hydraulic fracture propagation in deep reservoirs are discussed. Finally, the proposed model is extended to simulate the simultaneous propagation of multiple hydraulic fractures.

Model Verification
Since no analytical solution of elasto-plastic hydraulic fracture propagation can be referred to, the correctness and effectiveness of the presented model and algorithm are verified from the following two aspects: (1) fracture propagation under the elastic deformation-KGD model, and (2) fracture propagation under plastic deformation-ABAQUS.

Fracture Propagation under Elastic Deformation
Neglecting the plastic deformation and matrix flow, the proposed model is simplified to the classical KGD problem [6,7]. The numerical results are now compared with the existing analytical solutions under elastic deformation. The input parameters are set as follows: E = 25 GPa, ν = 0.2, t c = 0.5 MPa, G C = 50 N/m, Q inject = 0.0005 m 3 /s, and µ = 0.1 Pa·s. Figure 5 shows the comparisons of the fracture half-length and the fracture width at wellbore between the numerical solutions and analytical solutions [45]. It can be found that the numerical results agree well with the analytical solutions.

Fracture Propagation under Plastic Deformation
Another case is carried out to simulate hydraulic fracture propagation in the deep reservoir (100 m × 100 m). As shown in Figure 6, only half of the domain is simulated based on the symmetry of the problem, and the hydraulic fracture propagates along the direction of the maximum horizontal principal stress. To verify the efficiency of the proposed method for modeling fracture propagation under plastic deformation, obtained solutions by the proposed EDFM-XFEM model are compared with the results of the finite element model proposed by Wang [31] and solved by ABAQUS, a commercial software. The grids used for the simulations are also shown in Figure 6

Fracture Propagation under Plastic Deformation
Another case is carried out to simulate hydraulic fracture propagation in the deep reservoir (100 m × 100 m). As shown in Figure 6, only half of the domain is simulated based on the symmetry of the problem, and the hydraulic fracture propagates along the direction of the maximum horizontal principal stress. To verify the efficiency of the proposed method for modeling fracture propagation under plastic deformation, obtained solutions by the proposed EDFM-XFEM model are compared with the results of the finite element model proposed by Wang [31] and solved by ABAQUS, a commercial software. The grids used for the simulations are also shown in Figure 6 The distributions of displacement and stress solved by the two models are compared, as shown in Figures 7 and 8. It can be found that the results based on the two different models show good agreement, which demonstrates the accuracy of the proposed EDFM-XFEM model for modeling fracture propagation under plastic deformation. The difference between the EDFM-XFEM model and the finite element model is that the latter must use the conforming mesh, while the former can use the non-matching mesh. Therefore, for the finite element model, the hydraulic fracture can only extend along the grid boundary and only straight fracture propagation can be simulated. The EDFM-XFEM model can simulate planar or non-planar hydraulic fracture propagation in structured grids without remeshing.
The distributions of displacement and stress solved by the two models are compared, as shown in Figures 7 and 8. It can be found that the results based on the two different models show good agreement, which demonstrates the accuracy of the proposed EDFM-XFEM model for modeling fracture propagation under plastic deformation. The difference between the EDFM-XFEM model and the finite element model is that the latter must use the conforming mesh, while the former can use the non-matching mesh. Therefore, for the finite element model, the hydraulic fracture can only extend along the grid boundary and only straight fracture propagation can be simulated. The EDFM-XFEM model can simulate planar or non-planar hydraulic fracture propagation in structured grids without remeshing.
The distributions of displacement and stress solved by the two models are compared, as shown in Figures 7 and 8. It can be found that the results based on the two different models show good agreement, which demonstrates the accuracy of the proposed EDFM-XFEM model for modeling fracture propagation under plastic deformation. The difference between the EDFM-XFEM model and the finite element model is that the latter must use the conforming mesh, while the former can use the non-matching mesh. Therefore, for the finite element model, the hydraulic fracture can only extend along the grid boundary and only straight fracture propagation can be simulated. The EDFM-XFEM model can simulate planar or non-planar hydraulic fracture propagation in structured grids without remeshing.  The distributions of displacement and stress solved by the two models are compared, as shown in Figures 7 and 8. It can be found that the results based on the two different models show good agreement, which demonstrates the accuracy of the proposed EDFM-XFEM model for modeling fracture propagation under plastic deformation. The difference between the EDFM-XFEM model and the finite element model is that the latter must use the conforming mesh, while the former can use the non-matching mesh. Therefore, for the finite element model, the hydraulic fracture can only extend along the grid boundary and only straight fracture propagation can be simulated. The EDFM-XFEM model can simulate planar or non-planar hydraulic fracture propagation in structured grids without remeshing.

Results and Discussion
Next, a series of numerical cases are performed to investigate the influences of rock plasticity, internal friction angle, dilatancy angle and permeability on hydraulic fracture propagation. In addition, the proposed model is extended to simulate multiple hydraulic fracture propagations.

Evolution of Plastic Zone
As shown in Figure 8, it can be found that the area with the large stress change appears near the fracture tip, so the stress concentration occurs at the fracture tip. Figure  9 displays the evolution of the plastic zone during elasto-plastic hydraulic fracture propagation for the simulation case in Section 4.1.2. It can be observed that plastic zones appear around the fracture tip owing to the stress concentration. Besides, with the extension of the hydraulic fracture, unloading occurs behind the advancing tip and new

Results and Discussion
Next, a series of numerical cases are performed to investigate the influences of rock plasticity, internal friction angle, dilatancy angle and permeability on hydraulic fracture propagation. In addition, the proposed model is extended to simulate multiple hydraulic fracture propagations.

Evolution of Plastic Zone
As shown in Figure 8, it can be found that the area with the large stress change appears near the fracture tip, so the stress concentration occurs at the fracture tip. Figure 9 displays the evolution of the plastic zone during elasto-plastic hydraulic fracture propagation for the simulation case in Section 4.1.2. It can be observed that plastic zones appear around the fracture tip owing to the stress concentration. Besides, with the extension of the hydraulic fracture, unloading occurs behind the advancing tip and new plastic zones are formed around the new tip, so the plastic zone moves forward. Meanwhile, plastic strain accumulates around the fracture walls as illustrated in Figure 10.

Effect of Plastic Deformation
To analyze the effects of plastic deformation on the hydraulic fracturing process, the simulation case in Section 4.1.2 is performed again without considering rock plasticity. Figure 11 shows the comparisons of the fracture width profiles, fracture width at wellbore, fracture half-length and net pressure at wellbore. Figure 11a displays the comparison of the fracture width profiles for the same fracture length. The result shows that the opening of the elasto-plastic hydraulic fractures is larger than that of the elastic hydraulic fractures. Besides, as can be seen, the initial fracture length is memorized in the fracture width profile for the elasto-plastic hydraulic fracture, since no plastic deformation is associated with the initial fracture (as shown in Figure 10). The comparisons of the fracture width at wellbore and the fracture half-length are presented in Figure 11b and Figure 11c, which indicate that the elasto-plastic hydraulic fractures are wider and shorter than the elastic hydraulic fractures at the same injection time. Figure 11d demonstrates that the fluid pressure required for the propagation of elasto-plastic fractures is higher than that of elastic fractures. As rock plasticity results in higher energy dissipation in deep reservoirs, larger pressure is required for fracture extension, which leads to hydraulic fractures with larger openings. Meanwhile, corresponding to the mass conservation law, a larger opening means a smaller length, so the extension speed of elasto-plastic fractures is lower than that of elastic fractures. Therefore, to obtain the desired fracture length, larger injection volumes and higher injection pressures are needed for hydraulic fracturing in deep reservoirs.

Effect of Plastic Deformation
To analyze the effects of plastic deformation on the hydraulic fracturing process, the simulation case in Section 4.1.2 is performed again without considering rock plasticity. Figure 11 shows the comparisons of the fracture width profiles, fracture width at wellbore, fracture half-length and net pressure at wellbore. Figure 11a displays the comparison of the fracture width profiles for the same fracture length. The result shows that the opening of the elasto-plastic hydraulic fractures is larger than that of the elastic hydraulic fractures. Besides, as can be seen, the initial fracture length is memorized in the fracture width profile for the elasto-plastic hydraulic fracture, since no plastic deformation is associated with the initial fracture (as shown in Figure 10). The comparisons of the fracture width at wellbore and the fracture half-length are presented in Figure 11b,c, which indicate that the elasto-plastic hydraulic fractures are wider and shorter than the elastic hydraulic fractures at the same injection time. Figure 11d demonstrates that the fluid pressure required for the propagation of elasto-plastic fractures is higher than that of elastic fractures. As rock plasticity results in higher energy dissipation in deep reservoirs, larger pressure is required for fracture extension, which leads to hydraulic fractures with larger openings. Meanwhile, corresponding to the mass conservation law, a larger opening means a smaller length, so the extension speed of elasto-plastic fractures is lower than that of elastic fractures. Therefore, to obtain the desired fracture length, larger injection volumes and higher injection pressures are needed for hydraulic fracturing in deep reservoirs.
between the EDFM-XFEM model and the finite element model is that the latter must use the conforming mesh, while the former can use the non-matching mesh. Therefore, for the finite element model, the hydraulic fracture can only extend along the grid boundary and only straight fracture propagation can be simulated. The EDFM-XFEM model can simulate planar or non-planar hydraulic fracture propagation in structured grids without remeshing.
The distributions of displacement and stress solved by the two models are compared, as shown in Figures 7 and 8. It can be found that the results based on the two different models show good agreement, which demonstrates the accuracy of the proposed EDFM-XFEM model for modeling fracture propagation under plastic deformation. The difference between the EDFM-XFEM model and the finite element model is that the latter must use the conforming mesh, while the former can use the non-matching mesh. Therefore, for the finite element model, the hydraulic fracture can only extend along the grid boundary and only straight fracture propagation can be simulated. The EDFM-XFEM model can simulate planar or non-planar hydraulic fracture propagation in structured grids without remeshing. Further, the sensitivities of the two plastic parameters involved in the elasto-plastic tangent stiffness tensor Dep, internal friction angle φ, and dilatancy angle ψ are studied. The numerical simulation case in Section 4.1.2 is conducted again with the parameters remaining unchanged, except the internal friction angle which is set as 30°, 35°and 40°. Figure 12 displays the numerical results for the fracture width, fracture length and net pressure with different values of φ. As can be seen, with the increase in φ, the fracture width and injection fluid pressure decrease while the fracture half-length increases, which means the effect of rock plasticity on hydraulic fracturing is weakened as the internal friction angle increases. Then, to study the effects of the dilatancy angle ψ, the internal friction angle is fixed at 30° while the values of ψ are set as 10°, 15° and 20°. The numerical results are provided in Figure 13. The results show that the fracture width and injection fluid pressure decrease while the fracture half-length increases with the decrease in the dilatancy angle, which means the effect of the dilatancy angle on hydraulic fracturing is opposite to that of the internal friction angle. Compared with the latter, the former has little influence on hydraulic fracturing. Further, the sensitivities of the two plastic parameters involved in the elasto-plastic tangent stiffness tensor D ep , internal friction angle ϕ, and dilatancy angle ψ are studied. The numerical simulation case in Section 4.1.2 is conducted again with the parameters remaining unchanged, except the internal friction angle which is set as 30 • , 35 • and 40 • . Figure 12 displays the numerical results for the fracture width, fracture length and net pressure with different values of ϕ. As can be seen, with the increase in ϕ, the fracture width and injection fluid pressure decrease while the fracture half-length increases, which means the effect of rock plasticity on hydraulic fracturing is weakened as the internal friction angle increases. Then, to study the effects of the dilatancy angle ψ, the internal friction angle is fixed at 30 • while the values of ψ are set as 10 • , 15 • and 20 • . The numerical results are provided in Figure 13. The results show that the fracture width and injection fluid pressure decrease while the fracture half-length increases with the decrease in the dilatancy angle, which means the effect of the dilatancy angle on hydraulic fracturing is opposite to that of the internal friction angle. Compared with the latter, the former has little influence on hydraulic fracturing. between the EDFM-XFEM model and the finite element model is that the latter must use the conforming mesh, while the former can use the non-matching mesh. Therefore, for the finite element model, the hydraulic fracture can only extend along the grid boundary and only straight fracture propagation can be simulated. The EDFM-XFEM model can simulate planar or non-planar hydraulic fracture propagation in structured grids without remeshing.

Effect of Permeability
In this part, we investigate the influences of permeability on the hydraulic fracturing process. The numerical simulation case in Section 4.1.2 is carried out again with the input parameters unchanged, except the matrix permeability km, whose value varies from 5 mD to 20 mD. Figure 14 displays the distribution of pore pressure in the matrix with different

Effect of Permeability
In this part, we investigate the influences of permeability on the hydraulic fracturing process. The numerical simulation case in Section 4.1.2 is carried out again with the input parameters unchanged, except the matrix permeability k m , whose value varies from 5 mD to 20 mD. Figure 14 displays the distribution of pore pressure in the matrix with different permeabilities. It can be found that when the matrix permeability is low, a low pressure zone exists in front of the fracture tip. As the value of the matrix permeability increases, the pressure waves propagate outward in an approximate round shape and the low pressure zone gradually disappears. Figure 15 provides the variation in the fracture width at wellbore and the fracture half-length obtained with different values of matrix permeability. One can see that at the same injection time, as the value of matrix permeability increases, the fracture width and fracture half-length decrease. This is because, with the increase in matrix permeability, more fracturing fluid leaks into the matrix, which leads to a shorter and narrower hydraulic fracture.

Extension to Multiple Hydraulic Fractures Propagation
In the above cases, only a single straight fracture is simulated. To investigate the effect of rock plasticity on the propagation path of non-planar fractures, two hydraulic fracture propagations in poro-elastic media and poro-elasto-plastic media are performed in this part. The fracture spacing is 5 m. The other model settings and parameters remain unchanged, except the horizontal stress difference which is reduced to 5 MPa. Figure 16a,b show the comparisons of the fracture widths. It is obvious that when considering plastic deformation, hydraulic fractures with a larger opening and smaller length are obtained,

Extension to Multiple Hydraulic Fractures Propagation
In the above cases, only a single straight fracture is simulated. To investigate the effect of rock plasticity on the propagation path of non-planar fractures, two hydraulic fracture propagations in poro-elastic media and poro-elasto-plastic media are performed in this part. The fracture spacing is 5 m. The other model settings and parameters remain unchanged, except the horizontal stress difference which is reduced to 5 MPa. Figure 16a,b show the comparisons of the fracture widths. It is obvious that when considering plastic deformation, hydraulic fractures with a larger opening and smaller length are obtained, which is consistent with the conclusion of single hydraulic fracture propagation in Section 4.2.2. The comparison of fracture trajectories is shown in Figure 16c. It is interesting to find that the deflection of fracture direction in poro-elastic meida is smaller than that in poro-elastoplastic media, meaning that the stress-shadowing effect is enhanced when considering the plastic deformation. A possible explanation is that a larger fracture width leads to larger induced stress, which causes fractures to apply larger normal stresses to each other, thus the stress-shadowing effect is enhanced.

Conclusions
In this paper, a hybrid numerical model has been proposed to simulate elasto-plastic hydraulic fracture propagation in deep reservoirs. Based on the proposed model, the poro-elasto-plastic effects during hydraulic fracturing in the deep reservoir are investigated, as follows: 1. During hydraulic fracturing in deep reservoirs, rock plasticity results in higher energy dissipation in deep reservoirs, larger pressure is required for fracture extension, which leads to hydraulic fractures with larger openings. Therefore, larger injection volume and higher injection pressure are needed to obtain the desired fracture length; 2. As the internal friction angle increases from 30° to 40°, the fracture width and injection pressure decrease while the fracture propagation velocity increases, which means the effect of plastic deformation is weakened. The effects of the dilatancy angle on hydraulic fracturing are opposite to that of the internal friction angle; 3. As the matrix permeability increases from 5 mD to 20 mD, more fracturing fluid leaks into the matrix, leading to a shorter and narrower hydraulic fracture; 4. For two hydraulic fracture propagations in deep reservoirs, plastic deformation can enhance the stress-shadowing effect.

Conclusions
In this paper, a hybrid numerical model has been proposed to simulate elasto-plastic hydraulic fracture propagation in deep reservoirs. Based on the proposed model, the poroelasto-plastic effects during hydraulic fracturing in the deep reservoir are investigated, as follows:

1.
During hydraulic fracturing in deep reservoirs, rock plasticity results in higher energy dissipation in deep reservoirs, larger pressure is required for fracture extension, which leads to hydraulic fractures with larger openings. Therefore, larger injection volume and higher injection pressure are needed to obtain the desired fracture length; 2.
As the internal friction angle increases from 30 • to 40 • , the fracture width and injection pressure decrease while the fracture propagation velocity increases, which means the effect of plastic deformation is weakened. The effects of the dilatancy angle on hydraulic fracturing are opposite to that of the internal friction angle; 3.
As the matrix permeability increases from 5 mD to 20 mD, more fracturing fluid leaks into the matrix, leading to a shorter and narrower hydraulic fracture; 4.
For two hydraulic fracture propagations in deep reservoirs, plastic deformation can enhance the stress-shadowing effect.

Data Availability Statement:
The data presented in this paper are available on request from the corresponding author.