A Semi-Analytical Approach to Model Drilling Fluid Leakage Into Fractured Formation

Loss of circulation while drilling is a challenging problem that may interrupt operations, reduce efficiency, and may contaminate the subsurface. When a drilled borehole intercepts conductive faults or fractures, lost circulation manifests as a partial or total escape of drilling, workover, or cementing fluids, into the surrounding rock formations. Loss control materials (LCM) are often used in the mitigation process. Understanding the fracture effective hydraulic properties and fluid leakage behavior is crucial to mitigate this problem. Analytical modeling of fluid flow in fractures is a tool that can be quickly deployed to assess lost circulation and perform diagnostics, including leakage rate decline, effective fracture conductivity, and selection of the LCM. Such models should be applicable to Newtonian and non-Newtonian yield-stress fluids, where the fluid rheology is a nonlinear function of fluid flow and shear stress. In this work, a new semi-analytical solution is developed to model the flow of non-Newtonian drilling fluid in a fractured medium. The solution model is applicable for various fluid types exhibiting yield-power-law (Herschel-Bulkley). We use high-resolution finite-element simulations based on the Cauchy equation to verify our solutions. We also generate type-curves and compare them to others in the literature. We demonstrate the applicability of the proposed model for two field cases encountering lost circulations. To address the subsurface uncertainty, we combine the developed solutions with Monte-Carlo and generate probabilistic predictions. The solution method can estimate the range of fracture conductivity, parametrized by the fracture hydraulic aperture, and time-dependent fluid loss rate that can predict the cumulative volume of lost fluid. The proposed approach is accurate and efficient enough to support decision-making for real-time drilling operations.


Introduction
Naturally fractured formation is often prone to severe mud loss during drilling operations. As a result, problems may emanate because of the severity of lost fluid such as kicks, wellbore instability, formation damage, and sometimes freshwater aquifer contamination (Seyedmohammadi 2017). One procedure to encounter and mitigate this problem is to add lost circulation material (LCM) to the drilling fluid. LCM is commonly used in drilling applications to reduce and stop lost circulation (Attong, Singh, and Teixeira 1995;Ali, Kalloo, and Singh 1997;Olsen et al. 2019;Knudsen et al. 2015). The fluid characteristics of LCM, such as density and viscosity should be carefully selected based on the formation hydraulic properties such as the conductivity of the thief layers and fractures, among other factors (Luzardo et al. 2015). Due to the timescale of the problem, there is a need to develop accurate and efficient modeling tools applicable to realtime drilling operations to perform diagnostics and predictions. Field observations suggest that, during mud solutions for a particular fluid-type case. For general cases, we verify our solutions with full-physics numerical simulations. In the third section, we introduce new type-curves with the corresponding dimensionless groups. In section four, before the conclusion, we demonstrate the applicability of the proposed approach for two field cases.

Physical & Mathematical Model
The general governing equation used to describe non-Newtonian fluid dynamics is given by Cauchy equation (Irgens 2014;Cioranescu, Girault, and Rajagopal 2016), such that,  The shaded brown area shows the mud radial invasion, is wellbore radius, w is fracture aperture, is mud-loss volume.
Assuming steady-state conditions and neglecting gravity with low inertial effect in comparison to other forces, Eq. (1) becomes, (2) In 1D radial system, the above equation simplifies to (Panton 1984), Where ( , ) zr  is the radial shear stress component, perpendicular to the z -direction.
On the other hand, the shear stress component  is described by the Herschel-Bulkley fluid model (Hemphil, Pilehvari, and Campos 1993), that is, In Eq. (4), the parameter 0  is the yield shear stress, which determines the fluidity state, as explained later.
The flow index n is a positive number that reflects fluid rheological behavior. For instance, the fluid exhibits shear-thinning behavior when 1 n  , and shear-thickening when 1 n  . Typical values for flow behavioral index in drilling fluid ranges from about 0.3-1.0 (Kelessidis et al. 2006). The other parameters correspond to the consistency multiplier m , and the derivative of the radial velocity r v in the z-direction, reflecting the shear rate. Note that the Herschel-Bulkley model in Eq. (4) can also be used to describe Newtonian fluids when considering 0 0  = and 1 n = .

Solution Method
Our system of equations is given by the Cauchy equation (3) and the Herschel-Bulkley fluid model Eq.(4). The first equation combines two external forces, pressure force, and shear force, applied to the fluid. The second equation describes the fluid rheological behavior as a function of the two external forces. The proposed solution is based on the assumption that the drilling fluid viscosity is higher than the in-situ water viscosity in the fracture, with no significant mixing between the two fluids. Under these conditions, a pistonlike displacement is considered at the mud-water interfaces, which is a reasonable assumption (Razavi et al. 2017). Therefore, the fluid pressure and the shear force are approximated only within the mud invaded zone. The fracture is assumed infinite acting with constant average hydraulic aperture, and no-slip boundary conditions (B.C.) are set at the fracture walls.
The concept of the proposed solution method is illustrated in Fig. 2. As the fluid propagates deep into the fracture, radial flow velocity decreases, and shear stress diminishes. Therefore, shear-thinning is highest in the vicinity of the wellbore and gradually reduces with radial distance, exhibiting a shear thickening behavior. Furthermore, driven by the vertical variations of flow velocity and shear stress, fluid layers along the fracture aperture develop and introduce self-friction with the highest intensity at the fracture wall and reduces linearly towards the fracture centerline, as illustrated in Fig. 2 (4)). This region is denoted by the plug flow region, while the rest is the free flow region. As the fluid propagates further within the fracture, the plug region expands toward the fracture walls, eventually reaching a total stall of the flow, as illustrated in Fig. 3a. This condition represents the ultimate steady-state, where the pressure drop between the wellbore and the mud front becomes too small to overcome the yield stress 0  (see Fig. 3b).
In Eq. (5), plug z is the vertical extension of the plug region.
In the plug region, 0 r dv dz = , therefore, writing Eq. (6) for each region separately, one gets, The total volumetric flow rate total Q , is written as the sum of the rate within the plug region plug Q , and within the free region free Q , that is, On the other hand, the flux can be expressed in terms of the surface integral for each region by: 2 1 1 1 2 2 2 n tot n nn l n a dp Q dp dp dp d r wn r dr n w w w m n n n dr dr Be rearranging Eq. (10), it can be simplified to an ODE with a quadratic form, as follows, 2 0 n tl n ota BD r Q dp dp dr dr A where, the quantities ,, AB and D are nonlinear functions of n , m , w , as discussed in Appendix A.
By solving for the quadratic equation (11) Finally, by integrating Eq. (12) and recalling the total flux, we reach the final ODE system to solve: The ODE system in Eq. (13) is nonlinear and cannot be solved analytically, except for a particular case when 1 n = (see Appendix B). Therefore, a numerical ODE solver is used to solve Eq. (13) starting with initial conditions ( ) In Appendix B, we show the derivation of a particular case of the proposed solution when 1 n = . Fig. 4 compares the solutions of the mud invasion front (dimensionless) versus a dimensionless time obtain by our model and the one proposed by Liétard et al. (2002), which shows identical match. The solutions are presented for various values between 0.002 and 0.04 of the dimensionless parameter  , which is defined by, This test case shows that a particular solution ( 1 n = ) of the proposed semi-analytical model converges to the analytical solution of Liétard et al. (2002) for Bingham plastic fluids.

Model Verification: General Solution.
We use simulations to verify our semi-analytical solutions for general Herschel-Bulkley fluids with different flow indexes, n . Simulations were performed using a finite-element method within COMSOL-Multiphysics (Littmarck and Saeidi 1986) to solve the flow of a non-Newtonian fluid within two parallel circular plates, mimicking the geometry of the radial fracture, shown in Fig. 1. We emphasize that the numerical solution is based on the Navier-Stokes equations, whereas the semi-analytical solution is based on the more general Cauchy equation of motion corresponding to the Herschel-Bulkley fluid. To address the discrepancy between the two solution methods, we use the approach by Papanastasiou (1987) , which links the Navier-Stokes equation with the Cauchy equation by manipulating the Herschel-Bulkley fluid as, Where, the effective viscosity as a function of shear rate, 9 p m , we found that p m =100 to 500 provided reasonable accuracy and robustness. High-resolution simulations were also needed to reduce instabilities and numerical artifacts. The accuracy of this numerical solution has been discussed previously (Albattat and Hoteit 2019).
In the simulation model, the fracture is assumed to be initially saturated in water at a constant pressure. The size of the fracture in the simulation domain is selected to be large enough to mimic infinite boundary conditions. The simulation solution includes the pressure solution in the mud and water zones, and the mud-front propagation (Albattat and Hoteit 2019). Fig. 5 compares the semi-analytical solutions and the numerical solutions obtained for two cases with 1 n = , and 0.7 n = , reflecting a yield stress shear-thinning fluid. The solutions describe the mud invasion distance versus time. The plateau section in the curves reflect the maximum invasion distance, that is when the mud-front stalls. Both solution methods are in good agreement.

Type-Curves
Type-curves are often generated with dimensionless groups to provide quick interpretations and diagnostics of time-dependent trends. This technique is commonly used in well testing (Lee 1982). In this work, we adopt the following dimensionless variables: the dimensionless time,  is a dimensionless parameter reflecting the fluid rheological and invasion behavior. Additional analysis of the effect of  is discussed below. The governing system of equations given in Eq. (13) is then rewritten in terms of the dimensionless variables, where the same solution method is applied. Additional details on the derivation of the dimensionless form of the semi-analytical solution are shown in Appendix C. Fig. 6 shows a set of type-curves generated for different flow indexes = 1.0, 0.8, 0.6, 0.4, and different values of  . The type-curves show cumulative mud loss increase versus time.
At far enough distance from the wellbore, the plugging effect becomes more significant, and the fluid eventually halts. We compare our derived type-curves with the one proposed by (Majidl et al. 2010 In this work, however, we kept a second nonlinear order in the equations. As a result, it was not possible to generate a full analytical solution, but rather a semi-analytical solution, as detailed in Appendix A.
We found that the accuracy of this simplified solution by Majidi et al. depends on the selected problem parameters.

Method Demonstration
The proposed semi-analytical method is applied to four wells from two different fields. In Field case 1, the drilling mud was Bingham plastic fluid with known fluid properties. A deterministic approach is used to match the data with the semi-analytical model. In Field case 2, a more complex Herschel-Bulkley fluid was used. In this case, a probabilistic approach is adopted to account for various uncertainties in the mud and subsurface formation properties. Alterations in the mud properties could occur for various reasons such as mud/water in-situ mixing, transient thermal effects, among others.

Field Case 1.
In the first case, a drilling Bingham plastic fluid was used for the drilling of two wells: Machar18 and Machar20, in the Machar field in the North Sea (Liétard et al. 2002). Lost circulations within naturally fractured formation were encountered in both wells, as shown in Fig. 9. The fracture apertures were estimated to be around 0.42 and 0.64 mm for Machar 18 and Machar 20, respectively, based on a simplified method proposed by (Huang, Griffiths, and Wong 2011). The limitation of Huang's et al. method is in the assumption that the total mud-loss volume must be known, which is available in this case. In our proposed approach, however, the total mud-loss volume can be predicted by fitting the transient mud-loss behavior versus dimensionless time.
We note that the well data, given in Fig. 9, as provided by Liétard et al. 2002, are plotted versus different dimensionless groups than the ones used here. For consistency, we replotted these field data using the proposed dimensionless variables given in Eq. (17), as shown in Fig. 10. The semi-analytical solutions were then generated to replicate the transient mud-flow, which produced an excellent match to the trends, as depicted in Fig. 10. In this case, the behavioral flow index is unity (n=1) for Bingham plastic fluid model, 19.5 lb/100ft 2 of yield value, and 35 cp of plastic viscosity. The dimensionless parameters are α=0.00215 and 0.0006436, respectively. Fracture apertures of Machar18 and Machar20 wells were found to be =0.425 mm and 0.616 mm.

Field Case 2.
The second field case corresponds to two wells in the Gulf of Mexico (Majidi et al. 2008;Majidl et al. 2010). The mud-loss volumes (gallons per minute) were reported versus time in the two zones for a limited period before the mud-loss stops, as shown in Fig. 11. The pressure drop, which is the difference between the injection pressure and initial formation pressure, was reported to be within the range of 700-800 psi, and mud properties are n=0.94, m=0.08 lbf.s/100ft 2 , and τ 0 =8.4 lbf/100ft 2 at surface conditions. Because of the nature of uncertainties in the drilling fluid and subsurface properties, we apply a probabilistic modeling approach using Monte-Carlo simulations. This probabilistic approach is needed to account for various uncertainties in the fluid property alterations related to subsurface temperature conditions and fluid mixing (Babu 1998;Rommetveit and Bjorkevoll 1997). We perform uncertainty analysis by combining our semi-analytical solutions with Monte-Carlo simulations. The whole process is computationally efficient and could be performed within seconds. We vary six input parameters, including the flow index, yield stress, fracture aperture, consistency factor, and pressure drop. Fig. 12 (right) corresponds to the parameters in Zone 2. These distributions are used as inputs for the semianalytical model, which is driven by the Monte-Carlo simulation process. The history-match corresponds to the time-dependent cumulative mud-loss volume, where uncertainties are taken into consideration. Fig. 13 (left) shows the solution matching band of the field data using the semi-analytical solution combined with Monte-Carlo simulations. The matching p10, p50, and p90 percentiles are also shown. Fig. 13 (right) uses a violin plot to highlight the distribution behavior of the solution versus time. Fig. 14 shows similar analyses done to match the data in Zone 2.
For Well 1, the obtained probabilistic predictions for the fracture aperture were 0.57, 0.79, and 1.01 mm for p10, p50, and p90, respectively. The predicted total mud-loss volumes were 968, 2591, and 5331 bbl for p10, p50, and p90, respectively. The probabilistic distributions of the matching parameters are shown in Fig. 15 for Well 1, and in Fig. 16 for Well 2. b) Distribution of total mud loss uncertainty. Fig. 16-Histograms showing the predicted probabilistic distribution of the fracture aperture (a) and the total mud-loss volume (b) for Zone 2.

Discussion
The objective of the proposed model is to develop a numerical diagnostic tool to help predict the dynamic behavior of mud leakage into fractured formation. The tool can be used to perform quick sensitivity analyses and what-if scenarios to optimize LCM selection. The model is based on simple computations and runs efficiently in a spreadsheet. This feature makes it convenient to be coupled with Monte-Carlo simulations to address uncertainties, where thousands of simulations are typically needed. The input data for the model corresponds to the drilling fluid rheological properties, the injection pressure-drop, and the mud-flow rate. The modeling procedure consists of 1) converting and plotting the cumulative mud volumes versus a dimensionless time, and 2) matching the trend with the semi-analytical model or the type-curves by tuning the fracture hydraulic aperture. As a result, the matching trend provides an estimate of the fracture hydraulic aperture and a prediction of the mud-loss behavior, including the flow halt time and the total ultimate volume of mud leakage.
It should be noted that the proposed model is based on several simplifications, and therefore, its applicability is conditional. The main simplifications include the assumption of mud-flow in a horizontal fracture. This assumption was adopted to neglect the buoyancy effect, driven by the difference of densities of the drilling mud and the in-situ water. Consequently, the model is not expected to be applicable when the gravity effect is significant such as in the case of mud with a specific gravity different from one, flowing in inclined fractures. The model also neglects the effect of water displacement, downstream of the mud. Therefore, the pressure of the water phase is assumed constant, equal to the initial pressure of the formation. This assumption is valid when the mud viscosity is significantly higher than the water viscosity. Note that an analogy to this assumption is in Richards equations for water flow in unsaturated porous media, where the air pressure is assumed constant. Other simplifications in our model include ignoring the mud/water in-situ mixing and thermal effects. Furthermore, the modeling procedure may produce nonunique solutions. This problem is common in all curve-fitting approaches, such as pressure transient analysis (PTA) and decline curve analysis (DCA), which mostly occur if the trend of field data is not well established. Therefore, uncertainty analysis should always be accounted for, and different methods should be used to confirm the results.

Conclusions
Lost circulation during drilling operations is a common problem that requires immediate intervention to circumvent fluid loss. Diagnostic tools, based on simplified input data such as fluid properties, pressure, and rate trends, can be quickly deployed to quantify uncertainties related to the fluid leakage into the subsurface formation and to perform predictions. The main contributions of this work are the following: • A new semi-analytical approach is presented to model the leakage behavior of general Herschel-Bulkley fluids into a horizontal infinite-acting fracture, mimicking the effect of a fractured formation. • The analytical approach is based on simplified assumptions such as horizontal fractures with uniform aperture. • The proposed solution is a generalization of other analytical solutions developed for Bingham plastic fluids. The model is applicable to different types of non-Newtonian fluids, including yield stress shear-thinning and shear-thickening fluids. • The model was verified for general cases using high-resolution finite element simulations.
• The modeling approach can predict the trend of mud leakage in a system with horizontal fractures as a function of time. It can predict the effective hydraulic aperture of the fracture, the ultimate total mud-loss volume, and the expected duration before the leakage stalls, if conditions allow. • We introduced dimensionless groups and generated type-curves, which can provide quick diagnostics about the leakage behavior from matching the type-curve trends without a need for simulations. • We demonstrated the applicability of the model for four wells from two different fields. A numerical procedure was described to couple the model with Monte-Carlo simulations to perform predictions under uncertainties. • The proposed semi-analytical model is based on simple calculations that can be performed efficiently. This model is useful as a numerical diagnostic tool to perform quick predictions and to help optimize LCM selection by performing what-if scenarios.

Appendix A: Derivation from Cauchy Equation to transient fluid loss
The Cauchy equation of Motion is given by: At the steady-state without inertial and gravity effects, we get, In a radial system 1D, the above equation simplifies to (see (Panton 1984)): Eq.(20) relates the pressure to shear stress. Another equation to relate shear stress to velocity is given by Herschel-Bulkley model, such that, Eq.(23) is valid along the fracture aperture in the z-direction, where two regions, according to the velocity profile, can be distinguished: plug (non-deformed) region and free (deformed) region. We can see from Herschel-Bulkley equation model that the plug region corresponds to zero derivative of the velocity in the z-direction, that is, The velocity is defined at the three B.C. as, Substituting, we find the total flux, 1/ 1 1/ 0 0 0 1/ 4 2 2 1 1 2 n n total n dp n Q dp dp dp dr n dr dr d Simplifying the above expresses to get, To avoid obtaining complex number because of the negative dp dr inside the power term, we rewrite the equation as, The last term is approximated by using the second-order Taylor expansion, that is, 2 0 0 0 4 11 2 2 1 1 1 2 2 2 n tota n l n n dp Q r dp dp dp dr dr dr d n r w n n w w w m n n n  which can be solved as (Choi, Yun, and Choi 2018;Vidunas 2008;Ismail and Pitman 2000), We substitute back these values into Eq.(46) to get, We now have an analytical solution of pressure as a function of radial distance r . The Eq. (48), (45) and (44) are used when the total flow rate entering the fracture is known. Otherwise, if the total flow rate is unknown, we can use the total flow rate as, Implementing into Eq.(55) to generate a dimensionless differential form, This is an ordinary differential equation that can be solved analytically by using Inverse Function Theorem (IFT) in such a way that any function ( ) fx is both differentiable and invertible. Assuming In the above equation, the polylogarithmic function Bingham Plastic fluid model. Furthermore, an analytical solution was provided as well (Liétard et al. 2002(Liétard et al. , 2002.