The fractured carbonate rock reservoir is widespread in Sichuan Basin, and the characteristics of different areas are different. The development of natural fractures is varying degrees, lost circulations occur frequently, and the formation heterogeneity is strong, which causes the formation not sufficiently stimulated by acidizing. It may affect the effectiveness of reservoir stimulation. To advance the whole stimulation effect of the heterogeneous fractured carbonate reservoir, a new solution for determining the invasion radius during drilling in a fracture network reservoir is presented, which is based on solute transport and convection-diffusion equations. It can predict the invasion radius caused by mud loss and determines the range of mud loss invasion, which clarify the scope and degree of reservoir damage. The formation of skin factors polluted by mud loss was calculated. The experiments verified that the acidizing technology can remove reservoir damage and reduce the polluted formation of skin factors. The opening pressure of the nature fracture closed is calculated which can control the acidizing area. It is confirmed how many fractures in the carbonate reservoir can be opened under the wellhead pressure limit, which meets the construction conditions of acidizing fractured reservoirs. The framework of network-fracture deep acidizing technology was established, which can efficiently break through the detrimental zone caused by lost circulation, break down the natural fracture network, and decrease the formation of the skin. The restart pressure of natural fractures was calculated, and the design parameters such as pump pressure and displacement were optimized to quantify the scope of reservoir stimulation and the scale of acid fluid. The technique of network-fracture deep acidizing was applied for well A, the formation of skin after acidizing can be reduced to -4, and the testing production of well A was 58.87×104m3/d. The technique of network-fracture deep acidizing can quantify the acid scale and sweep area in acid fracturing design, which develops the fracturing efficiency and improves the fracturing engineering.

Fractured reservoirs are widely distributed in Sichuan Basin. In recent years, exploration practice has confirmed that fractured reservoirs have great exploration and development potential [14]. The fractured reservoir has low porosity and permeability, strong heterogeneity, and locally developed natural fractures and vugs which cause drilling fluid leakage during drilling and completion; therefore, formation damage is serious. Given the development of natural fractures in this kind of reservoir, the phenomenon of lost circulation is common in the drilling process, and the reservoir is polluted and blocked to varying degrees, so the main goal should be to remove plugging acidification and dredge the natural fracture cave system [58].

Due to the development of natural fractures in the reservoir, drilling fluid leakage is common, and the drilling fluid invades the reservoir with the natural fracture system to form a large gel leakage damage zone; the conventional plugging removal acidification process is difficult to make the acid fluid break through the leakage pollution zone and cannot fully remove the reservoir pollution blockage [915]. Therefore, it is necessary to predict the depth of the damage zone formed by the invasion of drilling fluid into the reservoir, which can be achieved by measuring the changes in formation properties. Researchers have mainly studied the invasion model of drilling fluid filtrate to predict the invasion depth range and verified the model through field data fitting. At the same time, it can also be used to measure and monitor the properties of near-wellbore zones and reservoirs [1618]. The drilling fluid intruded into the near-wellbore formation is displaced and mixed with the formation fluid, resulting in the formation of a certain depth of reservoir damage zone near the near-wellbore formation [1922]. Civan proposed an improved migration equation of two-phase fluid in deformable porous media taking into account the diffusion, cake formation, and filtrate invasion of compressible/incompressible cake model with/without solid invasion [2327]. Ramakrishnan and Wilkinson established a radial invasion model of water-based drilling fluid filtrate, which can determine the saturation of oil and filtrate in the formation water phase [28]. Pordel Shahri et al. simulated the convection-diffusion movement of drilling fluid filtrate in fractures with limited conductivity and the convection-diffusion movement facing the foundation block through the fracture wall [29]. The depth of damage zone formed by drilling fluid invading the reservoir is closely related to factors such as formation properties (porosity, permeability, gas saturation, and electricity), engineering parameters (drilling differential pressure and soaking time), and drilling fluid performance (density, water loss, and its system). For field applications, although too many parameters will make the prediction depth more accurate, the calculation is more complex, and the model calculation parameters may not be all available. Therefore, it is necessary to develop a new, fast, and accurate method to predict the depth of mud leakage invasion.

Another key factor is to determine how many closed natural fractures are fractured and treated by acidizing. It is also important to judge the fracture deformation behavior caused by acidizing. Hydromechanical coupling models were established based on the phase field method to investigate the effect of vugs on hydraulic fracture propagation, which were solved by finite element approximation and iterative method for the hydraulic fracturing problem [3034]. Numerous kinds of literature have been carried out to analyze the influence of acid fracturing on formation rocks. Fu et al. presented an explicitly integrated, fully coupled discrete finite element approach for the simulation of hydraulic fracturing in arbitrary fracture networks, and their model was verified against the K–G–D closed-form solution for the propagation of a single hydraulic fracture and validated against laboratory testing results on the interaction between a propagating hydraulic fracture and an existing fracture [35]. Ghommem et al. validated a 3D core scale predictive model for carbonate acidizing and showed a good capability to qualitatively capture the experimentally observed wormhole morphology and quantitatively predict the pore volume to break through [36]. Ma et al. developed a three-dimensional thermal-hydro-chemical (T-H-C) coupled unified pipe-network method (UPM) to simulate the process of acidizing treatment in fractured carbonate rock, which considered thermal effects and incorporated multiphysical governing equations for explicitly represented fractures in the current method [37]. Guo et al. further carried out a series of hydraulic fracturing experiments and found that etching wormholes induced by the rock-acid reaction would change the pore structure of carbonate rocks, and gelled acid enables to increase both the width and the roughness of main fracture surfaces [38]. Gou et al. conducted true triaxle fracturing experiments on carbonate rocks and observed that viscoelastic surfactant acid made the hydraulic fracture more conducive to connecting the nature fractures [39]. Jin et al. established and verified a damage constitutive model considering the compression hardening process of the acid-corroded sandstone under uniaxial loading, and proposed a theoretical model for evaluating the brittleness index (BI) of the acid-corroded sandstone based on energy evolution theory and damage constitutive relation [36]. Cao et al. proposed a method to dissect the internal structure of fault-karst reservoirs that is helpful for y targeted acid fracturing, which improved the prediction accuracy and the recovery efficiency of fault-karst reservoirs [41]. Xing et al. studied the uniaxial compressive strength (UCS) and uniaxial tensile strength (UTS) by five types of rock samples and obtained the fracture extension degree of each lithology via acid dissolution experiments [42].

By establishing the leakage and migration model of drilling fluid in the fracture network, the leakage and invasion depth of drilling fluid is quantified, the range of reservoir damage zone caused by well leakage is defined, and the natural fracture opening pressure is calculated. Based on not forming new fracturing fractures, the natural fracture system is fully blocked and dredged, the reservoir skin coefficient is reduced, and the network-fracture depth acidification technology based on leakage and invasion depth and natural fracture opening pressure is formed. Through field application, the scope of reservoir leakage damage zone of fractured formation in well A is predicted, and the opening pressure of natural fractures in the reservoir is calculated, to optimize the acidizing design parameters of the well, improve the acidizing efficiency, and achieve good acidizing transformation effect on site.

2.1. Reservoir Physical Properties

Taking the fractured reservoir in Sichuan China, Western Sichuan as an example, its average depth is 6700 m, the average temperature is 150°C, the average formation pressure is 125 MPa, the average porosity is 0.75%, and the average permeability is 0.95 mD (Figure 1). The observation of the core and microthin section reveals that the reservoir space of the fractured formation is mostly fractures and corrosion holes associated with fractures. The relationship between core porosity and permeability is poor, and most samples have the characteristics of low porosity and relatively high permeability, indicating that fractures play a very important role in improving reservoir porosity and permeability (Figure 2).

2.2. Rock Mechanical Properties

Triaxial compression tests were carried out on rocks of the Maokou formation in Western Sichuan, Sichuan Basin. The experimental results show that the average compressive strength of the Maokou group rocks is 445.65 MPa and the average Young’s modulus is 3.60×104MPa, with an average Poisson’s ratio of 0.23 (Table 1). It can be seen from Figure 2 that there are visible natural fractures in the Maokou formation rock sample, which will inevitably lead to the reduction of the compressive strength of the rock sample.

The results of the differential strain test show that the maximum horizontal principal stress of the Maokou group rocks is 163.51 MPa, the minimum horizontal principal stress is 141.21 MPa, and the vertical principal stress is 190.66 MPa (Table 2). The horizontal two-dimensional stress difference is 22.3 MPa, and the horizontal stress difference coefficient is 0.16. Combined with the three-dimensional principal stress distribution law, it is considered that the formation of vertical fractures is the main form of fracturing in the Maokou formation. Combined with the paleomagnetic experiment and differential strain experiment (Table 2), it is concluded that the direction of horizontal maximum principal stress is NE125°.

3.1. Technical Connotation

Based on the geological and engineering characteristics of the reservoir, firstly, predict the depth of drilling fluid leakage and invasion, combined with the prediction model of reservoir skin factor, clarify the scope and degree of reservoir damage, and quantify the scope of reconstruction and the scale of acid fluid. Calculate the restart pressure of natural fractures, and optimize the design parameters such as pump pressure and displacement. Finally, improve the efficiency of reservoir stimulation of fractured reservoirs. The technique design process is shown in Figure 3.

3.2. Prediction of Drilling Fluid Leakage Damage Zone

For the fractured reservoir with natural fractures, the drilling fluid leakage is large and the reservoir damage is serious (Figure 4). The goal of reservoir stimulation is to remove the reservoir pollution and blockage and dredge the natural fracture cave system. Therefore, the key to reservoir stimulation is to break through the damage zone formed by the leakage of drilling fluid along the natural fracture system, unblock the natural fracture system, and reduce the reservoir skin coefficient to the greatest extent [43, 44].

3.2.1. Mathematical Model Based on Drilling Fluid Solute Transport and Convection-Diffusion

Due to the leakage and invasion of drilling fluid, drilling fluid with different concentrations is distributed in the natural fractures of the formation. The fluid at the open end of the borehole wall fracture is pure drilling fluid, which is regarded as it is a dimensionless concentration of 1. The dimensionless concentration of drilling fluid is consumed with the radial extension of the fracture in the process of fracture leakage and flow and percolates and diffuses on the fracture wall. At the same time, it is also diluted by the formation fluid, which gradually reduces its concentration, until the concentration at the maximum invasion depth of drilling fluid is close to 0; as shown in Figure 5, the yellow curve represents the dimensionless concentration change curve of drilling fluid in the fracture. Therefore, the drilling fluid is regarded as a pure fluid different from the formation fluid. Based on the solute convection-diffusion theory, the dimensionless concentration distribution of drilling fluid leakage flow in the fracture network is simulated and studied, and the dimensionless concentration distribution profile of drilling fluid in the fracture is used to calculate the maximum depth of drilling fluid leakage invasion. Based on the fracture network geological model and the description of drilling fluid leakage and invasion, the leakage and migration model of drilling fluid in the fracture network is established.

The leakage migration of drilling fluid in fractures includes two modes: one is the convection-diffusion migration of drilling fluid leakage flow in the fracture extension direction with high conductivity, and the other is the filtration diffusion of drilling fluid filtrate in the foundation block system perpendicular to the fracture wall. It is assumed that the drilling fluid flows stably and uniformly in the radial fracture with a smooth fracture wall, and the corresponding control equation of the model is
(1)nfCft=Df2Cfr2+1rCfrvCfrq+qw,(2)rwrre,t>0,
where nf is the fracture volume fraction in the fracture network, Cf is the drilling fluid concentration in the fracture network (M/L3), Df is the diffusion coefficient of drilling fluid concentration in radial fracture (L2/T), r is the penetration depth of drilling fluid leakage (L), w is the average half-width of fracture (L), q is drilling fluid leakage rate (L3/T), q is drilling fluid leakage rate drilling fluid convection exchange flux between fracture and base block (L3/T), v is steady flow velocity of drilling fluid in radial fractures (L/T), re is the maximum penetration depth of drilling fluid leakage (L), and t is leakage time (T).
The calculation formula of fracture porosity in the fracture network model is
(3)nf=h1H1+h2H2,
where h1 and H1 are radial fracture width and fracture spacing (L) and h2 and H2 are plane fracture width and fracture spacing (L).
The velocity of the linear flow of drilling fluid at the open end of the fracture at the borehole wall is
(4)vrw=q2πrwh,
where vrw is the leakage rate of drilling fluid at the wellbore (L/T), q is the leakage rate of drilling fluid (L3/T), rw is the e radius of the wellbore (L), r is the radial invasion radius of drilling fluid (L), and h is the effective thickness of leakage layer (L).
The concentration diffusion coefficient D of drilling fluid consists of convection-diffusion coefficient De and molecular diffusion coefficient D, as shown in the following formula:
(5)D=De+D,
where De is the convection-diffusion coefficient of drilling fluid (L2/T) and D is the molecular diffusion coefficient (L2/T).
The hydrodynamic diffusion coefficient of drilling fluid in a fracture network is related to the leakage rate and diffusion of drilling fluid, and the convection-diffusion coefficient De can be expressed as the empirical relationship of radial velocity v of drilling fluid:
(6)De=αv,
where α is the diffusion factor of drilling fluid (L).
The drilling fluid diffusion exchange flux q between the fracture and the rock base block can be solved by Fick’s first law:
(7)q=DmnmϕmCfzi,nm=1nf,Dm=τD,τ1=1.61ϕm1.15,
where Dm is the diffusion coefficient of drilling fluid in rock base block (L2/T), nm is the volume fraction of the rock base block in the fracture network, ϕm is matrix porosity of rock base block, τ is the tortuosity of natural fracture, and Cf/zi is the concentration gradient of drilling fluid on the contact surface between fracture and base block.
t=0, drilling fluid filter cake is not formed at the well wall, and at this time, the concentration gradient of drilling fluid is the maximum. The filter cake of drilling fluid gradually forms with the leakage time, and the concentration gradient of drilling fluid gradually decreases to 0, when the drilling fluid filter cake is completely formed. This pressure gradient depends on the concentration of drilling fluid in the fracture, which can be calculated by the following formula:
(8)Cfzi=ACFCfeBCFt,
where ACF and BCF are the first constant and the second constant of the drilling fluid concentration on the fracture base block contact surface, respectively, which can be determined by indoor experimental fluid and rock. At the initial moment, t=0, Cf/zi is the maximum value and equals to ACFCf.
The exchange flux q between the fracture and the rock base block caused by convection is
(9)q=VfmCf,vfm=KmnmϕmμPzi,
where Vfm is the flow rate of drilling fluid on the indirect contact surface between fracture and rock base block (L/T), μ is the viscosity of drilling fluid (M/L. T), and Km is the permeability of rock base block (L2).
The pressure gradient on the contact surface between fracture and rock base block is
(10)Pzi=APFPf¯PieBPFt,
where Pf¯ is the fluid pressure in the fracture (Pa), Pi is the original formation pressure (Pa), and APF and BPF are the first constant and the second constant of the drilling fluid pressure on the fracture base block contact surface, which can be determined by experimental fluid and rock.
Substituting equations (3)~(10) into equation (1), the final form of the differential equation of drilling fluid leakage and invasion is obtained:
(11)q2πrhCDfr+αq2πrh2CDfr2+1rCDfr+DmnmϕmwACFCDfeBCFtKmnmϕmμwAPFPf¯PieBPFtCDf=nfCDft.
The dimensionless concentration CDf of drilling fluid in the fracture network can be expressed as
(12)CDf=CfCo,
where Co is the concentration of drilling fluid in the wellbore (M/L3).

3.2.2. Model Boundary Condition

(1) Model Initial Condition. When the drilling fluid does not leak into the fracture network, the dimensionless concentration of drilling fluid in the fracture network is
(13)CDf=0rw<r<re,t=0.
(2) Model Boundary Condition. The dimensionless concentration of drilling fluid at the borehole wall is
(14)CDf=1r=rw,t>0.
The dimensionless concentration of drilling fluid at the maximum depth of drilling fluid leakage invasion is
(15)CDfx=0r=re,t>0.

3.2.3. Model Calculation

Using the block central grid system, the finite difference method is used to separate the difference of equation (11), the central difference of distance r, and the forward difference of time t, and the fully implicit discrete form is obtained:
(16)aiCDfi1t+1+biCDfit+1+ciCDfi+1t+1=CDfit,
where ai=q/2πrihΔt/2Δr+αΔt/Δr21/riΔt/2Δr, bi=1+2αq/2πrihΔt/Δr2DmnmϕmΔt/wACFeBCFtn+1+kmnmϕmΔt/μwAPFpf¯pieBPFtn+1, ci=q/2πrihΔt/2ΔrαΔt/Δr2+1/riΔt/2Δr.

Formula (16) is to solve the dimensionless concentration distribution of drilling fluid with the invasion depth of drilling fluid leakage in the process of drilling fluid leakage.

The difference of the initial condition equation (13) of the model is obtained:
(17)CDfi=1i=0,(18)CDfi=0i=1,2,,n.
The model boundary condition equation (14) is differentiated to obtain:
(19)CDf0=1,CDfn=CDfn1=0,t=0,1,2,,TMAX.
Combined with equation (17) and equation (18), a matrix equation group is formed to solve the distribution of dimensionless concentration of drilling fluid with the invasion depth of drilling fluid leakage in the process of drilling fluid leakage:
(20)b1c1a2b2c2an1bn1cn1anbnCDf1t+1CDf2t+1CDfn1t+1CDfnt+1=CDf1ta1CDf0t+1CDf2tCDfn1tCDfntcnCDfnt+1.

Using mathematical software MATLAB, we can calculate and solve the matrix equations of drilling fluid leakage invasion depth distribution, and then, the distribution of drilling fluid leakage and invasion depth in natural fractures can be determined.

3.3. Calculation of Restart Pressure of Natural Fracture

3.3.1. Mechanical Model of Fractured Formation with Natural Fracture

Firstly, two cases of natural fracture opening are discussed, and the model shown in Figure 6 needs to be established.

The W-T criterion for tensile fracture is equation (21), and the Mohr-Coulomb criterion for shear fracture is equation (22).
(21)pσn,(22)pσ0σμf+σn,
where p is the pore pressure near the wall of natural fracture (MPa), σn is normal stress acting on the natural fracture surface (MPa), σ is shear stress acting on the natural fracture surface (MPa), σ0 is the cohesion of rock in natural fracture (MPa), and μf is friction coefficient of the natural fracture surface dimensionless.
After the artificial joint intersects and connects with the natural joint, the liquid enters the natural joint, and its near wall face pressure p is
(23)px,t=σh+pnetx,t.
The net fracture pressure required for tensile fracture is
(24)pnetx,tσHσh21cos2θ.

The tensile fracture of natural fractures needs to meet ΔPPnetΔσ.

For shear failure, the net fracture pressure required for shear fracture is
(25)pnetx,t1μfσ0+σHσh2μfsin2θμfcos2θ,
where Pnet is net pressure in joint (MPa).

It can be seen that the shear failure of natural fractures needs to meet ΔPPnetΔσ+σ0/μf.

3.3.2. Mechanical Model of Fractured Formation without Natural Fracture

For reservoirs without natural fractures, the above two criteria are not applicable, and new models and mechanical criteria need to be established. The mechanical model of fracture formation in a rock mass with undeveloped fractures is shown in Figure 7.

Elliptical fractures are formed in the plane wireless domain, the long half axis is lf, the short half axis is w, axial compressive stress at infinity is σH, and compressive stress in the short axis direction is σh, and there is uniform pressure in the seam Pnet.

The boundary condition of the model:

When y=0,xlf,σy=Pnet,σxy=0, and x2+y2, we can take as σxσH,σyσh,σxy0.

According to the elastic failure criterion of the rock body, the net fracture pressure is obtained:
(26)pnet=σHσh+St.
where St is the tensile strength of rock (MPa).

To break the rock body and produce new joints, it is necessary to meet PnetmaxσHσh+St; that is, temporary plugging of natural fractures does not develop, and the reservoir forms new branch fractures in the rock mass need to meet ΔPPnetSum of maximum and minimum horizontal stress difference and rock tensile strength.

3.3.3. Critical Conditions for Closed Natural Fracture Opening

The stress of natural fractures depends on the total stress state of the formation and the orientation and dip angle of the fractures. Even if the natural fractures in the stratum are at the same depth, if they develop in different directions and dip angles, their normal effective stress is also different. For horizontal joints and vertical joints, they are only affected by vertical effective stress and horizontal effective stress, respectively. It is assumed that the principal stresses on the crack surface along the X, Y, and Z axes are σx, σy, and σz.

The approach angle between natural fracture and shaft or hydraulic fracture is θ, and the normal stress σn and shear stress τn on the natural fracture surface are obtained by
(27)σn=σH+σh2+σHσh2cos290θ,(28)τn=σHσh2sin290θ,
where θ is the approach angle between natural fracture and shaft or hydraulic fracture (°), σH is the maximum horizontal principal stress (MPa), and σh is the horizontal minimum principal stress (MPa).
Suppose that a small amount of acid liquid percolate into the natural fracture before the natural fracture opens. If the fluid pressure drop in the natural fracture is not considered, the fluid pressure in the natural fracture is equal to the bottom hole pressure. Assuming that the normal stress on the wall of the natural fracture is σn, when the fluid pressure in the natural fracture is greater than the normal stress σn acting on the natural fracture, the originally closed natural fracture will open. Therefore, the criterion of natural fracture opening is
(29)pf>popen=σn,
where pf is fluid pressure in the fracture (MPa) and popen is the natural fracture opening pressure (MPa).

The opening of natural fractures depends on the effective stress of natural fractures, which is not constant. When the fluid pressure in the natural fracture increases, the effective stress acting on the natural fracture surface decreases, and the fracture opening will increase. As the crack opening changes, the flow capacity of the fracture also changes, showing stress sensitivity.

To judge whether the closed natural fracture can open under a certain acid injection displacement, a simplified physical model is established as follows: inject liquid into the wellbore at a certain flow rate. When the injection volume is greater than the outflow volume, the wellbore internal pressure will increase due to fluid compressibility. When the bottom hole pressure exceeds the critical pressure of the natural fracture opening, it is judged that the natural fracture can open and meet the conditions of network-fracture acidizing construction.

Establish the material balance equation:
(30)ΔV=VinjVout,Pf=Pi+ΔP,
where ΔV is the volume change of fluid in unit time (m3), Vinj is the volume of injected fluid per unit time (m3) Vout is the volume of fluid flowing out per unit time (m3), Pi is the initial bottom hole pressure (MPa), and ΔP is the added extra pressure (MPa).
Assuming that the wellbore volume V is a certain value, according to the definition of fluid compressibility, ΔP can be expressed as
(31)ΔP=ΔVVwCw,
where Cw is the acid compressibility (MPa-1).
Assuming that the natural fracture is annular with the wellbore as the center, and the filtration of acid into the natural fracture and reservoir conforms to Darcy’s law, ignoring the improvement of reservoir permeability by acid in a short time, the amount of acid filtration into the reservoir per unit time is
(32)Vout=2πKhΔPlnrf/rw,
where K is the rock permeability of the perforated section (mD), H is the thickness of the perforated section (m), rf is the length of natural fracture (m), and rw is the borehole radius (m).
According to formulas (27)–(32), the critical displacement for opening natural fractures can be calculated. Considering the influence of the wellhead pressure limit, it is also necessary to determine the wellhead pressure under different acid injection displacements:
(33)pt=pf+pFph,
where pt is wellhead pressure (MPa), pF is the friction of acid in the wellbore (MPa), and ph is the liquid column pressure in the wellbore (MPa).

3.4. Simulation and Calculation of Reservoir Skin Factor

According to the calculation method of pollution damage skin coefficient Sd of the porous reservoir,
(34)Sd=KKd1lnrdrw.
The pollution radius of mud loss is calculated by the volume method:
(35)rd=Vlossπhϕ+rw2,
where rd is mud contamination radius (m), Kd is the average permeability of the contaminated zone (mD), Vloss is mud leakage (m3), and ϕ is the average porosity of the reservoir (%).
Based on Prats’ calculation method of skin coefficient of a hydraulic fracturing well, the effective action distance of acid fluid after network-fracture acidification and the change of natural fracture width is equivalent to the effective wellbore radius; the skin factor calculation model for predicting the acidizing effect of network fractures is established to guide the optimization of acid injection displacement and acid injection volume. After acidizing the fractured reservoir, the effective borehole radius is
(36)rew=refrwDa,
where ref is the effective action distance of acid (m), rew is the effective borehole radius (m), and rwD is the dimensionless equivalent borehole radius. Dimensionless parameters a can be defined as
(37)a=πkref2kfw=6πkrefw3,
where kf is the formation permeability (mD) and w is the fracture width (m).
When the acidification area is larger than the pollution radius, the acidification skin coefficient is
(38)Sacid=kkdlnrdrewlnrdrw.

3.5. Experimental Verification and Evaluation

The reservoir core is used for artificial fracture. During the acid displacement experiment, by increasing the confining pressure, it is ensured that the cut fracture remains closed during the displacement experiment. The experimental results show that under the condition of not opening the crack, the acid liquid does not penetrate the core, and only intense dissolution occurs locally at the end face of the core (Figure 7). On the other hand, by appropriately reducing the confining pressure, it is ensured that the split fracture remains open during the displacement experiment. The experimental results show that when the fracture is opened, the acid can easily break through the pollution zone and form a certain depth of etching grooves on the fracture wall (Figure 8). The experimental simulation results also show that for the reservoir with closed natural fractures, the acid solution mainly dissolves along the acid injection end face, and no wormhole is formed. For reservoirs with open natural fractures, acid first flows along the natural fractures and erodes the wall of the natural fractures (Figure 9), and wormholes expand at the tip of the natural fractures (Figure 10).

4.1. Prediction of Leakage Intrusion Depth in Well A

4.1.1. Model Parameter Assignment

The well section is 3330.4-3372.5 m. The gas reservoir in Section 1 is comprehensively interpreted, and the reservoir thickness is 35.9 m. The reservoir porosity is 1.1% ~6.9%, the average porosity is 2.5%, and the average permeability is 3.62 mD. Type I reservoir is 1.3 m, with an average porosity of 6.7%. Class II reservoir is 17.4 m, with an average porosity of 3.4%. Class III reservoir is 17.3 m, with an average porosity of 1.9%. The average number of fractures is 20/m, and the fracture width range is 0.17-0.41 mm. The basic parameters required by the model are shown in Table 3.

4.1.2. Numerical Simulation Results

There was one lost circulation and one gas invasion in the drilling of the oil test section (3330.0~3358.0 m), with a cumulative loss of 629.4 m3 and a fracture width of 1.7~4.1 mm according to logging interpretation. The reservoir geological parameters and field leakage data are substituted into the prediction model, and the prediction value of drilling fluid leakage invasion depth is obtained through a finite difference solution. According to the prediction results of numerical simulation (Figure 11), the invasion depth of drilling fluid increases with the lost circulation time. When the lost circulation occurs for 50 hours and the loss is 629 m3, the radius of the lost circulation damage zone of the drilling fluid is about 40 m.

4.2. Calculation of Natural Fracture Opening Net Pressure

According to the rock mechanics experiment in this well, the calculated cohesion is 24.1~92.4 mPa, with an average value of 69.3 mPa. The natural fractures with approach angle ≤60~80 are mainly tensile (Figure 12). According to the model calculation, the natural fracture opening PNET is 4.9~12.7 mPa (Figure 13), which provides a basis for the pump pressure in the stimulation construction.

4.3. Construction Pressure Calculation of Natural Fracture Opening

According to the formulas (27)–(33), the variation law of bottom hole pressure under different displacement conditions is calculated.

Calculate the construction friction, and calculate the pipe string concerning the friction: tubingΦ88.9mmpipewallthickness12.09mm×1800m+tubingΦ88.9mmpipewallthickness9.53mm×4670m+tubingΦ73mmpipewallthickness5.51mm×383m.

The formula for calculating the string friction is as follows:
(39)ΔPF=32ρq2fLπ2D5,
where ΔPF is the construction string friction (MPa), Ρ is the acid density (kg/m3), q is the pump displacement (m3/min), L is the well depth (m), and D is the tube diameter (m).

According to formula (39), take three points in the stable section of the curve to calculate the average friction factor F, and then, bring the calculated friction factor into formula (39) to calculate the friction and friction coefficient under different displacements. The calculation results are shown in Figure 14.

When the acid injection displacement is greater than 4.0m3/min, the conditions for natural fracture opening are met, and the displacement at this time is the minimum displacement required for natural fracture opening in formation. When the acid injection displacement reaches 5.0m3/min, the wellhead pressure is close to the pressure limit (125 MPa), and the displacement can no longer be increased. At this time, the bottom hole pressure is 160 MPa, and the displacement at this time is the maximum displacement that can be achieved during construction (Figure 15). Under deep conditions, the minimum fracture opening pressure is 143 MPa (0° inclination), 149.8 MPa (30° inclination), 155.4 MPa (60° inclination), and 165 MPa (90° inclination). Therefore, under this condition, 0°, 30°, and 60° natural fractures can open, while 90° natural fractures cannot open. It is confirmed that most of the fractures in the reservoir of this well can be opened under the condition of a wellhead pressure limit, which meets the construction conditions of acidizing fractured reservoirs.

4.4. Acidizing Design and Construction Optimization Analysis

4.4.1. Optimization of Acidizing Construction Parameters

According to the results of the prediction model of drilling fluid leakage invasion depth established above, it is concluded that the invasion depth of drilling fluid is approximately 30 m, so the acid corrosion fracture and acid fluid sweep range should be greater than 30 m in the acidizing design, FracPT software is used to simulate that when the scale of acid fluid is 320m3, and the corresponding length of acid corrosion fracture is 35 m, which can remove the leakage pollution of drilling fluid, as shown in Figures 16 and 17. Through the simulation of the reservoir skin coefficient, it is obtained that the reservoir skin coefficient after 320m3 gel acid acidification is -4 (reservoir depth is about 20 cm), which is shown in Figure 18. Combined with the prediction model of drilling fluid leakage and invasion depth, it shows that reservoir acidification in this scale can reduce the skin coefficient of the fractured reservoir.

4.4.2. Analysis of Acidizing Construction Effect

When the oil pressure of the acid solution in the well is blocked, the oil pressure of the acid solution in the well is gradually reduced, which means that the oil pressure of the acid solution in the well is blocked, and the oil leakage of the reservoir is gradually reduced in the acid injection stage, which means that the acid solution in the well is blocked, and the oil pressure of the reservoir is gradually reduced in the acid injection stage, as shown in Figure 19. The phenomenon of oil pressure drop in acidizing construction curve is analyzed, which shows that the acid solution relieves the reservoir pollution and blockage caused by lost drilling fluid and dredges the natural fracture cave system. After acidizing, the well obtained a high test production of 58.87×104m3/d. The field application proves that the prediction model of drilling fluid leakage and invasion can optimize the acidizing design of reservoir stimulation and improve the efficiency of reservoir stimulation.

For the fractured reservoir with serious leakage pollution, first, remove the pollution blockage near the well, reduce the reservoir skin coefficient, break through the leakage blockage zone, realize deep acidification, expand the transformation scope to increase the discharge area, and improve the efficiency of conventional plugging removal and acidification transformation.

Based on the drilling fluid leakage invasion model of the upper fracture network, the scope of the leakage damage zone is defined. Combined with the calculation of the restart pressure of natural fractures, the opening pressure of nature fracture closed is calculated which can control the acidizing area. It is confirmed how many fractures in the carbonate reservoir can be opened under the wellhead pressure limit, which meets the construction conditions of acidizing fractured reservoirs.

The formation of skin factors polluted by mud loss was calculated. The experiments verified that the acidizing technology can remove reservoir damage and reduce the polluted formation of skin factors.

The network-fracture depth acidification design parameters are further optimized, which provides a theoretical basis for field application.

The framework of network-fracture deep acidizing technology was established, which can efficiently break through the detrimental zone caused by lost circulation, break down the natural fracture network, and decrease the formation of the skin. The restart pressure of natural fractures was calculated, and the design parameters such as pump pressure and displacement were optimized to quantify the scope of reservoir stimulation and the scale of acid fluid. It has been applied to gas wells and achieved a remarkable acidizing effect.

Some or all data, models, and code generated or used during the study appear in the submitted article. All data, models, or code generated or used during the study are available from the corresponding author by request.

The authors declare no conflict of interest.

Mingwei Wang was responsible for the conceptualization, funding acquisition, project administration, resources, and software. Song Li was responsible for the data curation, methodology, formal analysis, methodology, and writing of the original draft and wrote, reviewed, and edited the manuscript. Wen Wu was responsible for the investigation, methodology, software, and visualization. Tao Li wrote, reviewed, and edited the manuscript. Gensheng Ni was responsible for the visualization and supervision. Yu Fu and Wen Zhou were responsible for the funding acquisition and wrote, reviewed, and edited the manuscript.

This work was supported by the National Natural Science Foundation of China: Study on Dynamic Characteristics of Methane/Carbon Dioxide in Shale Heterogeneous Reservoir under Multi-Field Coupling (41772150).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).