New Leakage Model for Naturally Fractured Formations Considering the Effects of Dual-System Hydro- Mechanical Coupling

Within the existing leakage model accounting for drilling mud loss in naturally-fractured formations, the leak-off velocity is assigned to a fixed value or described by the Cater model, which does not consider the influence of dual-system hydromechanical coupling effects between fracture-wall and fracture-inner systems. The dual-system between the formation and fracture is controlled by the flowing net pressure inside the fracture, which determines the dynamic width of the natural fracture and leak-off velocity. In this study, first, the leak-off velocity under the hydromechanical coupling of the fracture-wall system was obtained based on the coupled governing equations of the solid and liquid phases of the natural fracture-wall, as well as Darcy’s law. Second, the leakage-front invasion velocity, leakage rate, and leakage volume under the hydromechanical coupling of the fracture-inner system were clarified according to the geometric governing of the natural fracture morphology. Finally, the dual-system coupling leakage model was developed considering the continuous equation, while the numerical solution was obtained through a time-step deduction. Results show that at a given time, a greater formation permeability leads to a greater leakage rate and volume, with a smaller leakage front distance. The leakage rate increases with an increase in formation permeability, well bottom differential pressure, and initial width of the natural fracture, while it decreases with an increase in the fracture normal stiffness, yield stress, and plastic viscosity. The new leakage model and numerical method concerning time-step deduction are assessed by solving the issues of fully coupled fracture-wall and fracture-inner systems considering drilling fluid leak-off. The new model may be utilized to simulate various problems related to the invasion of drilling fluids into the fractures, including predicting the dynamic width of natural fracture and borehole ballooning/breathing phenomena.


Introduction
According to statistics, the economic loss caused by lost circulation in the global oil industry exceeds hundreds of millions of dollars every year [1], and the leakage plugging cost of fractured formations accounts for more than 90% of the total cost of lost circulation [2]. For example, the time-loss caused by lost circulation during the drilling of ultradeep wells in western Sichuan is 8%. By analyzing the drilling fluid leakage data and understanding drilling mud leakage, we can judge and evaluate the fracture development characteristics in the formation and provide the required parameters for subsequent leakage plugging operations [3][4][5][6][7][8]. When drilling mud loss occurs in naturally fractured formations, Sanfillippo et al. proposed an analytical model using Newtonian fluid to describe the leakage, neglecting the effect of the rheological properties of drilling mud [9]. Lietard et al. and Majidi et al. developed models for the radial flow of non-Newtonian fluid into a single natural fracture, neglecting the effect of drilling fluid leak-off and fracture deformability [10][11][12]. Other researchers have proposed an elastic equation to incorporate the impact of the fracture deformability, which relates the fracture width variation to the local pressure through the fracture normal stiffness [13][14][15][16][17][18][19]. To analyze fluid leak-off, Bychina conducted threshold processing for this leak-off process; that is, when the well-bottom differential pressure is greater than a fixed value, the fluid seeps into the formation at a certain leak-off velocity [20]. Razavi et al. described this with the Cater model [21]. Their research shows that for permeable fractured formations, considering the impact of the fluid, the leak-off effect is necessary for a realistic characterization of the leakage. However, the above model does not consider the porous elastic media characteristics of naturally permeable fractures and the hydromechanical coupling effect of the fracture-wall. When drilling mud loss occurs in a naturally fractured formation, the initial flowing net pressure field inside the fracture causes fracture deformation and leak-off, and reaction pressure is generated. Due to the flowing net pressure, the leak-off of the drilling fluid inside the fracture permeates into the formation region close to the fracture-wall, changing the formation pore pressure and redistributing the corresponding stress field. The time-dependent disturbed stress field arising from the coupled hydromechanical process changes the natural fracture width. The dynamic width of natural fractures is determined by changes in the deformed stress field of the fracture-wall and correspondingly influences the variation of the flowing net-pressure field inside the fracture. This is the joint effect of dual-system hydromechanical coupling between the fracture-wall and fracture-inner systems. In this work, it develops a new leakage model based on the poroelastic characteristics of the fracture-wall system and geometric governing equations of the fracture-inner system. The model considers the hydromechanical coupling effects between the aforementioned dual systems and provides a realistic characterization of the drilling mud loss in a naturally fractured formation. Besides, it also delineates the variations concerning the influence of multiple factors on the leakage state.

Hydromechanical Coupling of Fracture-Wall System
When drilling mud loss occurs in a naturally fractured formation, owing to drilling fluid leak-off, the hydromechanical coupling effect of the natural fracture-wall should be analyzed to understand the characteristics of the porous elastic medium and meet the governing equations of solid and fluid phases.
(1) Solid governing equations Stress equilibrium equation ignoring the body force.
where σ ij is the component of total stress tensor in MPa. Strain-displacement equation (positive strain indicates compression).
where ε ij is the component of the bulk strain tensor and u i,j is the component of the deformation tensor in m. Strain-stress-pressure equation where σ KK = σ xx + σ yy + σ zz and δ ij , E, υ, and p are the Kronecker symbol, Young's modulus in MPa, Poisson's ratio, and fluid pressure in MPa, respectively. Other constants include the bulk modulus K = E/3ð1 − 2vÞ in MPa, Biot coefficient α = 1 − K/K s , and K s is Particle volume modulus in MPa.
From the above equation, the stress-strain-pressure is where ε e is the volumetric strain defined by ε e = ε x + ε y + ε z , the shear modulus G is G = E/2ð1 + vÞ in MPa, and λ is the Lame coefficient defined by Considering the relation between the volumetric strain ε e and the solid displacement vector u Furthermore, substituting Eq. (6) into Eq. (5) yields (2) Fluid governing equations Mass conservation, Darcy's law, and the state equation of the pore fluid in the formation region next to the natural fracture-wall are porous medium and are expressed as follows: Mass conservation Darcy's law [22] is written as The state equation reads as In Eqs. (8) through (10), ρ is the fluid density in kg/m 3 , ϕ is porosity expressed as a percentage, v f ! is the fluid velocity vector in m/s, t is time in s, v s ! is solid velocity vector in m/s, μ is fluid viscosity in Pa·s, κ is the formation permeability in μm 2 , p is the fluid pressure in Pa, and c f is the fluid compressibility in MPa -1 .
As leak-off occurs only along the fracture-wall (along the Z-direction), it is described by radial coordinates; thus, where ξ = z/ ffiffiffiffiffiffiffiffiffi 4c Π t p , by substituting it into the above equation, one obtains (c is an undetermined coefficient) As shown in Figure 1, the natural fracture leak-off is divided into two areas, namely, the drilling fluid invasion area (area 1) and original formation area (area 2). The pressure distribution in the corresponding areas is represented by p 1 and p 2 , respectively. The fluid viscosities of the two areas are the drilling fluid viscosity μ l and the formation viscosity μ, respectively. The term sðtÞ indicates the contact surface of the two areas; where for area 1, c Ι = MðK + 4G/3Þκ/ðK + 4G/3 + α 2 MÞμ l , ϖ = z/ ffiffiffiffiffiffiffi ffi 4c Ι t p . Owing to the continuous pressure and mass distribution on the contact surface, the boundary conditions of the contact surface are given by Eqs. (19) and (20).

Geofluids
The pressure p 1 of area 1 and p 2 of area 2 are calculated by Eq. (18). Combining the boundary condition in Eq. (19), the pressure at the contact surface sðtÞ is where erf ðξÞ is the error function, erfc ðξÞ is the complementary error function, and erfc ðξÞ = 1 − erf ðξÞ. To ensure that Eq. (21) is applicable at all times, η can be set as an undetermined coefficient. The terms η and sðtÞ have the following relationship:

Geofluids
Based on Eq. (23), the dichotomy method is used to solve for the undetermined coefficient η, yielding c 1 . Therefore, the drilling mud leak-off velocity under the fracture-wall hydromechanical coupling is Letting c 1 = −2Δp/ ffiffi ς p and substituting it into Eq. (24) yields where ς is the harmonic coefficient of the leak-off velocity. As leak-off occurs in both the upper and lower sides of the fracture-wall, the total leak-off velocity equation is where u L is drilling fluid leak-off velocity in m/s.

Hydromechanical Coupling of Fracture-Inner System
When drilling mud loss occurs in a naturally fractured formation, this can be a penny-shaped natural fracture with an infinitely large radius [11,18,20,21] (see Figure 2(a)). The initial flowing net pressure inside the fracture p N ′ causes fracture deformation and leakage and correspondingly produces the reaction pressure p S ′. The superposition of pressures p N ′ and p S ′ further yields the flowing net pressure inside the fracture pðrÞ (see Figure 2(b)). The flowing net pressure inside the fracture pðrÞ acts on the natural fracture-wall, causing it to deform. The linear deformation rule yields Eq. (27) [16,17,19]. In Eq. (27), the natural fracture is assumed to have a nonzero fracture width (w 0 ) at zero overpressure, which implies that the natural fracture remains open at the formation of pore pressure [23][24][25].
where wðrÞ is the dynamic width of natural fracture in m, w 0 is the initial width of the fracture (natural fracture width at pðrÞ = 0) in m, and pðrÞ is the net pressure flowing inside the fracture in MPa. The term K Frac is the fracture normal stiffness in MPa/m. The flowing net pressure inside the fracture along the radial direction is [21] p r where Δp is the well bottom differential pressure in MPa.
CðtÞ is the dimensionless morphological index coefficient of the natural fracture, and r f ðtÞ is the distance to the leakage front in m. By analyzing Eq. (28), at the fracture opening (r = r w ), p ðrÞ = Δp is the well bottom differential pressure; at the leakage front distance (r = r f ), pðrÞ = 0 is the internal pressure in the original fracture. Both the morphologic index coefficient of natural fracture CðtÞ and the leakage front distance r f ðtÞ are time functions, which are hereinafter referred to as C and r f . Substituting Eq. (28) into Eq. (27), the radial distribution of the dynamic width of the natural fracture is When drilling mud loss occurs in a naturally fractured formation, the leakage front distance r f increases from r f n−1 to r f n (i.e., time ðn − 1ÞΔt → nΔt), and the increment of the drilling fluid loss volume ΔV n within Δt is divided into two parts: the leakage volume increment ΔV F n due to the deformation of the natural fracture and the increase in the distance to the leakage front, and the leak-off volume increment ΔV L n of the permeable fracture-wall (see Figure 3).
According to the above analysis, ΔV n is obtained based on mass conservation where ΔV F n is solved by the following equation where r f n is leakage front distance at t = nΔt in m. Based on the leak-off velocity Eq. (26), ΔV L n is expressed as where ΔA is the increment leak-off area within Δt, m 2 .

Geofluids
The geometric governing equation of the leakage invasion velocity at any location r (r ∈ ½r w , r f where r f is the leakage front distance) is written as At this point, the total leakage volume V F ðtÞ due to the deformation of the natural fracture and increase in the distance to the leakage front is Meanwhile, the leak-off volume at the leakage distance r is the total leak-off volume within t minus the leak-off volume before leakage distance r; V L ðr, tÞ is expressed as By substituting Eqs. (29), (34), and (35) into Eq. (33), the geometric governing equation of the leakage rate is where the expression of Q is Substituting the leak-off velocity u L ðr * , tÞ from Eq. (26) into Eq. (36) yields where τ is the time required to reach the leakage distance r * . Substituting Eq. (38) into Eq. (33) yields the expression of the geometric governing equation of leakage invasion velocity v m r, t

Leakage Coupling Model for Dual Systems
When leakage occurs in a naturally fractured formation, it includes two sets of coupling systems for this process, the fracture-wall system (coupling system 1) and fracture-inner system (coupling system 2). The initial flowing net pressure field inside the fracture causes the leak-off and deformation of the fracture. Meanwhile, reaction pressure is generated. Due to the flowing net pressure, the leaked drilling fluid inside the fracture permeates into the formation close to the fracture-wall and changes the formation pore pressure; this causes the corresponding stress field to redistribute. The time-dependent disturbed stress field arising from the coupled hydromechanical process changes the natural fracture width. The dynamic width of a natural fracture is determined by changes in the deformed stress field of the fracture-wall and correspondingly reacts to variations of the flowing net-pressure field inside the fracture. The link between the two coupled systems of the formation and fracture is the flowing net pressure inside the fracture, which determines the dynamic width of the natural fracture and leak-off velocity (see Figure 4).  Figure 4: Schematic diagram of dual-system coupling.

Geofluids
When drilling mud loss occurs in a naturally fractured formation, the pressure gradient inside the fracture of the Bingham fluid model is [11] dp Substituting the geometric governing equation of the leakage invasion velocity Eq. (39) into Eq. (40), integrating r from r w ⟶ r f , and combining the flowing net pressure inside the fracture pðrÞ Eq. (28), the well bottom differential pressure Δp is written as After integrating Eq. (41), the well-bottom differential pressure Δp is given by where χ Ι = ð4πκΔp/μ l ffiffiffiffiffi ffi ςc Ι p Þ and the expressions of The pressure distribution inside the fracture can be obtained from Eq. (45) as To solve Eqs. (45) and (46), the leakage time t is divided into nΔt, that is, if t = nΔt, then r f = r f n . The leakage front distance increases from r f n−1 ⟶ r f n within ðn − 1ÞΔt ⟶ nΔt(see Figure 5). In this way, the spatial distance increment Δr f n corresponding to the time increment Δt has the same leak-off velocity, which allows the equation to be solved in segments. The item Ð r r w Þdr * dr in Eq. (46) is processed by segments. At this point, the leakage front distance is r f , and the upper limit of the integral is r ∈ ½r w , r f . The integrated domain is processed in segments and expressed as ∑ k=m k=0 Ð r f k+1 r f k f ðrÞdr, where m ∈ ½0, n − 1 is an integer; then, this item can be expressed as   Figure 6: Schematic diagram for solving the dual-system coupling leakage model. As shown in Eq. (28), performing segment manipulations of the expression of flowing net pressure inside the fracture yields the following expression of t = nΔt, pðrÞ written as By Eq. (36), at a given time, we can obtain the leakage rate at any leakage distance. The leakage rate measured on the ground is the leakage rate at the fracture opening, that is, qðr w , tÞ. Substituting the expression of Δp L into Eq. (36) and calculating it by segments, we obtain qðr w , tÞ q r w , t n ð Þ= Q dr f n dt + 1 where , ς k is calculated by the flowing net pressure inside the fracture of each spatial distance Δr f n in any period of time Δt (set Δt = 0:001 s). At a given time, the total leak-off volume is V L ðr w , tÞ = ∑ j=n j=1 q L j Δt, calculates it by segments provides V L ðr w , t n Þ Á C+1 Cr f k+1 + r f k+1 + r f j À Á C 2 + 3C + 2 As seen from Eq. (30), the total leakage volume measured on the ground is the sum of the two parts of V L ðr w , t n Þ and V F ðt n Þ The pressure inside the fracture solved by Eq. (48) is p I ðrÞ; the pressure inside fracture solved by Eq. (50) is p II ðrÞ; the rationality of C is determined by judging whether they satisfy the solution accuracy ε P = jp I ðrÞ − p II ðrÞj/p I ðrÞ. The whole solution process is divided into two steps "initial value determination" and "recursive prediction" (see Figure 6). The specific steps are as follows: (1) Enter the basic parameters Δp, τ y , μ p , w 0 , K Frac , and basic formation parameters (k, μ, E, ν, K s , and K f ) to calculate the leak-off harmonic coefficient ς 9 Geofluids rate qðt n Þ, leakage volume inside fracture V F ðt n Þ, leak-off volume V L ðt n Þ, and total leakage volume V Total ðt n Þ Among them, steps (2)-(4) represent the "initial value determination" step, and steps (5)-(7) are the "recursive prediction" step.

Method Verification and Leakage Simulation.
To analyze the leakage state at different formation permeabilities, κ was set as 0, 1, 5, and 10 md. Other parameters are shown in Table 1. When κ = 0 md, the calculation results were compared with Bychina's [20] finite difference method to verify the correctness of the new model. As shown in Figures 7-9, when κ = 0 md, the calculation results of the proposed model are consistent with those of Bychina's finite difference method, which verifies the reliability of the model. At any given time, a greater formation permeability causes a greater leakage rate and volume, as well as         a smaller leakage front distance. Although there is little difference in the leakage rate and leakage volume under different formation permeabilities, they have different physical meanings. With an increase in formation permeability κ, the leakage front distance decreases, which means that the leakage volume V F ðtÞ decreases and the leak-off volume V L ðtÞ increases. In addition, as shown in Eq. (38), with the increase of r f , the effect of κ on the leakage rate is evident. Moreover, the lower formation permeability produces a higher leakage invasion speed, which increases the difficulty of solid plugging [26].

Analysis of Multiple Influence Factors.
To analyze the influence of multiple influence factors under different           (1) Fracture normal stiffness As shown in Figures 10-12, at any given time, a greater fracture normal stiffness causes a smaller leakage front distance, leakage rate, and leakage volume under the same formation permeability. Based on Table 2, it can be seen that the formation permeability and fracture normal stiffness have the same effect on the leakage front distance, while having the opposite effect on the leakage rate and leakage volume. In other words, the leakage front distance decreases with an increase in the formation permeability and fracture normal stiffness. However, the leakage rate and leakage volume increase with increases in formation permeability and decrease with increasing fracture normal stiffness.
(2) Rheological properties of drilling fluid As shown in Figures 13-15, at any given time, a greater normal yield stress causes a smaller leakage front distance, leakage rate, and leakage volume under the same formation permeability. Based on Table 2, it can be seen that the formation permeability and yield stress have the same effect on the leakage front distance while having the opposite effect on the leakage rate and leakage volume. In other words, the leakage front distance decreases with an increase in the formation permeability and yield stress. However, the leakage rate and leakage volume increase with increasing the formation permeability and decrease with increasing yield stress.
As shown in Figures 16-18, at any given time, a greater plastic viscosity causes a smaller leakage front distance, leakage rate, and leakage volume under the same formation permeability. Based on Table 2, it can be seen that the formation permeability and plastic viscosity have the same effect on the leakage front distance while having the opposite effect on the leakage rate and leakage volume. In other words, the leakage front distance decreases with an increase in the formation permeability and plastic viscosity. However, the leakage rate and leakage volume increase with increasing formation permeability and decrease with increasing plastic viscosity.

Geofluids
As shown in Figures 19-21, at any given time, a greater well bottom differential pressure causes a greater leakage front distance, leakage rate, and leakage volume under the same formation permeability. Based on Table 2, it can be seen that the formation permeability and well bottom differential pressure have opposite effects on the leakage front distance, while having the same effect on the leakage rate and leakage volume. In other words, the leakage front distance decreases with an increase in the formation permeability and increases with an increasing well bottom differential pressure. However, the leakage rate and leakage volume increase with both an increasing formation permeability and well bottom differential pressure.
(4) Initial width of natural fracture As shown in Figures 21-24, at any given time, a greater initial width of the natural fracture causes a greater leakage front distance, leakage rate, and leakage volume under the same formation permeability. Based on Table 2, it can be seen that the formation permeability and initial width of the natural fracture have opposite effects on the leakage front distance while having the same effect on the leakage rate and leakage volume. In other words, the leakage front distance decreases with an increase in the formation permeability and increases with an increasing initial width of the natural fracture. However, the leakage rate and leakage volume increase with an increase in both formation permeability and initial width of the natural fracture.
In summary, the leakage front distance decreases with an increase in the formation permeability, fracture normal stiffness, yield stress, and plastic viscosity and increases with an increasing well bottom differential pressure and initial width of natural fracture. The leakage rate and leakage volume increase with increases in formation permeability, well bottom differential pressure, and initial width of the natural fracture; they decrease with increasing fracture normal stiffness, yield stress, and plastic viscosity. As shown in Table 2, the combined impact of an arbitrary combination of two or three sets of influencing factors on leakage front distance, leakage rate, and leakage volume can be obtained by multiplying the single influencing factor. When the impacts have opposite effects, the final result depends on the degree of influence of each factor.