Injection-production optimization of fault-karst reservoir—considering high-speed non-Darcy effect

Fault-karst reservoirs are featured by complex geological characteristics, and accurate and fast simulation of such kind of reservoirs using traditional simulator and simulation methods is pretty hard. Herein, we tried to discrete the complex fault-karst structures into one-dimensional connected units connecting the well, fracture and cave based on reservoir static physical parameters and injection-production dynamics. Two characteristic parameters, conductivity and connected volume, are proposed to characterize the inter-well connectivity and material basis. Meanwhile, the high-speed non-Darcy seepage term is introduced into the material balance equations for well-fracture-cave connected units to describe the actual seepage characteristics within the fault-karst reservoirs, and to better simulate the oil/water production dynamics. Based on this method, a fracture reservoir model of 1 injection-3 production was established. The change of oil–water action law in different injection and extraction systems under two production regimes of fixed production rate and fixed pressure is analyzed. A case study was conducted on S fault zone, where the flow of oil and gas did not follow the Darcy seepage rule and with a β value of 103–104, the single well flow pressure and oil production were perfectly matched with the real data. In addition, connected units with more prominent high-speed non-Darcy features were found to have better connectivity, which might shed light on the more accurate prediction of inter-well connectivity. Moreover, an improved injection-production well pattern and was proposed based on connectivity prediction model and reservoir engineering method to solve the problems of insufficient natural energy supply and overhigh oil production rate in Block S. Furthermore, the injection/production rate as well as the timing and cycle of water injection was predicted and optimized so as to better guide to site operations.


Introduction
In conventional oil and gas reservoirs, the seepage velocities of oil and gas phase are small, and the seepage processes satisfy the Darcy flow rule. For a typical Darcy flow, only the effect of viscous force is considered and the pressure gradient is proportional to the seepage velocity (Wen-Juan et al. 2014). Fault-karst reservoirs are mainly composed of caves and fractures with high permeability and porosity and matrix with extremely low permeability and porosity, in which case. fractures are the dominant seepage channels, and karst caves of different sizes are the major reservoirs for oil and gas. For pores and fractures with abundant oil and gas, if the oil/gas production rates are high, the seepage process would deviate from the linear Darcy flow rule (Wen-Biao et al. 2020). Therefore, in this paper, Forchheimer equation (Jiu-Gu and You-Sheng 2011) instead of Darcy equation was applied to more precisely describe the high-speed non-Darcy flow in faultkarst reservoirs. Meanwhile, node encryption method was used to depict the special structures of fractures and caves and to better reflect the reservoir features of fault-karst reservoirs. Herein, the complex seepage process in faultkarst reservoirs is discretized into connected unit based one-dimensional flow, and an oil/water dynamics based 1 3 interwell connectivity model was established on the basis of material balance equations and Berkeley displacement equations for non-Darcy flow.
The strong heterogeneity of fault-karst such kind of reservoirs poses a great difficulty for oil and gas exploitation. However, most of the current numerical simulation methods for fault-karst reservoirs have neglected the fact that the seepage flow in caves was no longer Darcy flow. Therefore, a novel simulation method that considered the impacts of high-speed non-Darcy seepage flow was proposed to more precisely simulate the dynamics of faultkarst reservoirs. Block S is a typical fault-karst reservoirs with karst caves as the main reservoir and cracks as the main dominant channel in Xinjiang. In Block S, only a few fault-karst units are in their stable elastic production stage, and most wells have already entered into the declining stage. With new wells continuously being put into production, the current production is on the rise. The high production rate at the early production stage generally leads to the rapid decline of natural energy and formation pressure, and oil production are to be terminated within a short period. The natural decline rate in Block S has increased significantly since 2019. Considering the high gas-oil ratio, high saturation pressure and high formation collapse risk in Block S, it is of urgent necessary to optimize the production rate, the injection-production well pattern and other operation parameters to solve or mitigate the problems occur in or follow the late elastic production stage. In this article, we developed an inter-well connectivity model for Block S to clarify the connecting paths within the reservoirs based on well-reservoir relationship, remainly oil distribution law and internal reservoir structure. Meanwhile, reservoir engineering methods were also introduced to design the optimization scheme for water injection and oil production so as to perfect the injection-production well pattern to balance the injection and production rate, to expand the sweep efficiency, and to guide the stable production.

Model establishment
(1) Establishment of material balance equation and pressure solving Discrete the well-pore-fracture-cave structure into nodes, and apply the Forchhmer equation (Macini et al. 2011;Huang et al. 2020;Yuan et al. 2013) to simulate the high-speed non-Darcy flow in fractures and caves by introducing a high-speed non-Darcy flow coefficient β into the material balance equation. Seepage velocity and pressure gradient follow the quadratic nonlinear relationship: and are the viscosity and density of fluid, respectively, k v and v v are the Forchheimer permeability and the velocity, respectively, and is the non-Darcy coefficient.
The material balance equation of fracture-cave unit in the interwell connectivity model which was studied by ZHAO-Hui (Hui et al. 2017(Hui et al. , 2014 was improved as: Where A ij is the interface area between control body i and control body j, m 2 ; L ij is the distance between well i and well j, m; P t i and P t j are the average pressures in the oil drainage area at time t, MPa; Q t i is the source sink phase of node i, m 3 /s; Q t ci is the channeling flow from fracture to karst cave, m 3 /s; C ti is the comprehensive compression coefficient, MPa −1 ; Δt is the time step, d; V ti is the single well control volume, 10 4 m 3 .
Dynamic conductivity of fluids that demonstrating non-Darcy flow features in fracture-cave structures was calculated based on the above-mentioned derivation process. By solving Eq. (2) with the implicit difference method, the average pressures in the drainage area of each well at different time were obtained (Hui et al. 2017(Hui et al. , 2014, and the fluid flowing direction as well as the flow rate between nodes could be determined. Figure 1 is the depiction of fracturecave connected units fault-karst; T ij is the interwell conductivity between wells i and j, representing the interwell oil/ water volume flow rate under unit pressure difference; V ij is the effective pore volume of the interwell connected units. where P t wfi , P t wfj are the bottom-hole flow pressure of production well and water injection well, respectively, Q t i , Q t wj are the Production and injection of wells I and j at time t, respectively, J t Oi is the oil production per unit production pressure difference, J t wj is the injection volume under unit injection pressure difference, N p is the cumulative oil recovery. (2) Saturation tracking (Hui et al. 2014) The oil/water saturations between wells and within fractures and caves were tracked based on the Berkeley-Villewater Water-flooding Theory. Take node X and its upstream node X u as the research objects, then: X is the seepage location, X u is the upstream of seepage point X, is the porosity, g ′ w is the derivative of water cut, S w , S wu are the water saturation at X and X u .
Take the influences of different measures implemented in actual reservoir production process into account and the derivative of interwell water cut at time t was obtained: where Q t Vij , Q t ′ Vij are the dimensionless cumulative flows at time t and t ' .
(3) Optimization of injection and production parameters 1) Selection of conversion well Clarify the characteristics of interwell connectivity based on the connectivity prediction model, select connected units with high connectivity and choose well groups with wellacknowledged reservoir developmental situation and high remaining oil saturation. As for a single well, reservoirs with constant volume were preferred and no impacts would be imposed on the neighboring reservoirs after the conversion. Moreover, it is preferable that the conversion well was surrounded by the cave-type and non-fractured reservoirs. The field experience showed that the well trajectory located at the bottom of the oil layer and in contact with the middle-upper part of the reservoir should be selected for water injection to construct a perfect bottom water displacement scheme.
2) Optimization of liquid production at water-injection stage Calculate the critical yield using the Dupuit formula (Yun 2017;Han et al. 2015;Zhao-Jie et al. 2015): where q c is the critical liquid production rate, K is the average permeability, L is the length of horizontal section, B o is the oil volume coefficient, Z w is the water avoidance height, r w is the well radius, is the reservoir anisotropy ratio, w , o are the density of formation water and oil, respectively, h is the thickness of the layer.
The average permeability is: where K h and K v are the horizontal and vertical permeability. The dimensionless water avoidance height: where h D is the dimensionless water avoidance height. The anisotropy ratio is: 3) Injection timing With the continuous production of oil, gas, and water, the formation pressure continues to drop, the interactions between fractures and caves become more obvious, and the seepage capacity between the fractures and caves decreases. When the formation pressure was lower than the critical fluid pressure at crack closure, crack plastic closure occurs. Therefore, the moment when the formation pressure equals to the critical fluid pressure is the best injection timing.
According to the calculation method of critical fluid pressure (Gen-Hua et al. 2021;Xiangjun et al. 2005), the effective normal stress of the fracture plane could be expressed as: where ne , H and h are the effective positive stress, maximum and minimum stress in the horizontal direction, v is the vertical principal stress, P p is the critical pressure.
Assuming that the micro-convex was isotropic and there were no interaction forces, then for a two-dimensional plane, there was: where x , y and z are the effective stress in X, Y and Z, is the Poisson's ratio of rock.
Drucker-Prager yield failure criterion: where z = ne = s .Where I 1 is the first stress invariant, m is the material parameter, K is the yield stress of rock. If F(σ) = 0, Substituting Eq. (14) into Eq. (11), the critical fluid pressure that satisfies the plastic yield conditions was:

4) Optimization of water injection rate
Consider the influences of water injection rate on displacement, and set the pressure recovery degree and stabilized water supply radius as the optimization objectives (Wen-Xue and Yong 2019; Wei 2017; Shen-Wen et al. 2015) for faultkarst reservoirs to quantitatively optimize the water injection volume. When the injection volume is too small, the gravity differentiation would lead to serious vertical flow of injection water, resulting in invalid water injection and energy loss; while when the injection rate is too high, the fast water pushing speed might lead to insufficient water displacement and obvious water channeling, resulting in a low sweep efficiency. Therefore, a suitable water injection rate is necessary to ensure the stable movement of the oil-water interface and to improve the sweep efficiency.
In the isolated elastic production stage, the oil production per unit pressure drop was calculated according to the material balance principle (Ai-Min 2007): where ΔP d is the pressure drop value, ΔP u is the pressure rising value.
To achieve the goals of stabilizing the flow conductivity and the supply radius of a single well on the basis of the current production system, the water injection rate should be: Couple Eqs. (16) and (18): If the supply radius did not reduce, and the injected water can sweep over the interwell recoverable remaining oil, then: where N is the formation crude oil reserves.

5) Optimization of injection cycle:
Water injection method affects the pressure interference range and degree in a given injection and production system. In order to increase the water sweep efficiency and to effectively prevent rapid water channeling, it is necessary to set a reasonable water injection cycle and as for cyclic injection, the well shut-in stage is the perfect time for water spreading in dominant fractures channel and oil production.

Application of the conceptual model
A conceptual model of fault-karst reservoirs was established based on the INSIM method and the FROTSIM simulator. The number of grids is 80 × 80, and grid sizes are 15 × 15 × 20 m 3 ; Fig. 2 shows the distribution of permeability, porosity, as well as the conductivity and connected volume of the connected units. There were 4 caves of different sizes randomly distributed in the reservoir, which were connected by fractures. The average permeabilities of caves and fracture were 800 mD and 1800 mD, respectively, and the pore volume of cave accounted for 97.18% of the total reservoir pore volume. The viscosities of oil and water were 40 MPa·s and 1 MPa·s, respectively. There are four wells in the model, among which P1, P2 and P3 wells were on the large karst cave, while W1 well was on the center of fracture.
The simulation process was divided into two stages: elastic production and water flooding stage. The elastic production lasted for 180 days, and the liquid production rates of P1, P2, P3 and W1 wells were 10 m 3 /d, 10 m 3 /d, 5 m 3 /d, 5 m 3 /d, respectively. The subsequent water flooding stage lasted for 1800 days, and the liquid production rates were 10 m 3 /d for P1 and P2 wells were, and 40 m 3 /d of injection rate for P3 well was.
As shown in Fig. 3, the fitting effect of oil production rate and water cut of wells P1, P2 and P3 is ideal, which verifies the correctness of the prediction model in this paper. Among them, well P1, well P2 and well P3 happened water breakthrough on the 1200, 900 and 1500 days, respectively.
The influence of high-speed non-Darcy flow term on conductivity was analyzed based on the connectivity prediction model, and the differences of production indexes at different β values were also compared.
(1) As shown in Fig. 4, the conductivity values of the 30th, 90th, 150th, 210th and 270th day in the elastic production stage, the conductivity values of the 330th, 630th, 930th, 1230th and 1500th day in the water flooding stage were taken for analysis.
As shown in Fig. 4, connectivity change of each connective connected units with at different production time was compared: When β = 0 (Darcy flow), the conductivity did not change with the production time. When β > 0 (non-Darcy flow), the conductivity of each connected unit decreased with the production time, and with the increase in β, the conductivity decreased faster and the non-Darcy flow characteristics become obvious.
Injection-production well: Fig. 4a shows the dynamic change of harmonic conductivity between W1 and P1. P1 (production well) was located on a center of karst cave; W1 (injection well) was located on a center of crack. In the cave-fracture connection system, there are great differences in porosity, permeability and structure between caves and fractures. The oil/water interaction is extremely complex, and the seepage potential at its boundary does not meet the linear distribution. Therefore, with the continuous production, the change rate of conductivity decreases significantly and decreases in a parabolic manner.
In Karst Cave: Fig. 4b shows the dynamic change of harmonic conductivity between P2 and V12. P2 (production well) was located on a center of karst cave; V12 was located on a boundary of crack. P2 well and V12 belong to karst cave system. The difference of geological characteristics between them is small, and there is no drastic change of flow potential. Therefore, with the production, the conductivity basically shows a linear downward trend.
Isolated fracture-cave: Fig. 4c shows the dynamic change of harmonic conductivity between V35 and V40. V35 was located on a boundary of crack; V40 was located on a center of karst cave. The conductivity in the isolated (2) As shown in Fig. 5, the influences of β value on the production dynamics of production wells P1, P2, P3 and block were analyzed.
P1 When the bottom-hole flowing pressure was a constant, the liquid, oil production rate and water cut increase with the β value increased. With the continuous production, when the β value is small, the liquid production rate shows a upward trend. But when the β value reaches 300, the liquid production rate shows a downward trend. And the water cut is lower than the condition under β = 100.
P2 When the bottom-hole flowing pressure was a constant, the liquid, oil production rate and water cut decline with the β value increased. With the continuous production, when the β value is small, the liquid production rate shows a downward trend. But when the β value reaches 300, the liquid production rate shows an upward trend.
P3 Always maintain constant liquid production rate. When the β value increased, the time of water breakthrough was later, the oil production rate increased and the water cut decreased.

Fig. 4 Connectivity change of each connected unit at different production time
Based on the dynamic change law of single well and block, the high-speed non-Darcy effect mainly occurs in the early and late stage of production. With the enhancement of high-speed non-Darcy effect, there is a critical value of β, which changes the law of oil-water interaction. There are differences in the law of oil/water interaction between constant liquid production rate and constant bottom-hole flowing pressure. It can be seen from Fig. 5j: as the relative permeability of water phase decreases, the oil phase basically remains unchanged, the oil/water mobility ratio increases, the oil phase seepage capacity is relatively enhanced, and the water phase is more obviously affected by (3) Figure 6 shows the production pressure difference of wells P1, P2 and P3. In the fixed liquid production rate stage, when the β value increased, the production differential pressure of the three wells increases in a fluctuating manner. Therefore, under the constant liquid rate system, when obtained the same production rate, the production pressure difference required to consider the high-speed non-Darcy effect is greater (Fang et al. 2016;Chuan-Zhi et al. 2011;Jing and Shi-Qing 2009;Ming-Yue et al. 2001). When P1 and P2 wells are converted to constant pressure condition, the production pressure difference decreases and increases in a fluctuating manner, respectively, with the β value increased.

Connectivity prediction
S fault zone is a fault-dominate karst-cave carbonate reservoir with normal temperature and pressure. The main oil/gas storage (Li-Xin et al. 2021) spaces in S zone are fractures and karst caves. Until December 2020, the cumulative oil production in S zone was 196.1655 × 10 4 t. At present, S zone is now in its stable production stage with low water cut. In order to simulate the zone more accurately, the proposed connectivity model was used to match the field production index and quantitatively predict the influences of high-speed non-Darcy free flow in the S-fault zone.
As shown in Fig. 7, the controllable reserves per well was divided and the control volume of each single well was obtained by historical matching. Table 1 shows that the reserve fitting results of each single well were good. Introduce the fitted parameters, including connected volume, average permeability and well spacing into Eq. 20 and the initial value of interwell conductivity was clear. Combine the automatic history fitting method and consider the high-speed non-Darcy flow features in the karst caves, the conductivity was obtained.
S fault zone is now still in its stable production stage and only a small amount of water was injected, the water content is low and the production is high. Herein, bottom-hole flow pressure was used as an important dynamic index, and a high-speed non-Darcy seepage term was introduced into the model to conduct better sensitivity analysis and to match the bottom-hole flow pressure; therefore, the influences of interwell Non-Darcy flow were quantified. Figure 8a illustrates that at Darcy seepage condition (β = 0), there were obvious errors between calculated and actual bottom-hole flow pressure values. When β = 10 3 -10 4 and β = 10 2 -10 3 , the fitting results of S-6 and S-2 bottom-hole flow pressures were perfect, respectively. Figure 8c is the sensitivity analysis and fitting results of cumulative oil production in block S. The fitting result was ideal when β was 10 3 . Meanwhile, since the zone is currently in a low water-cut stable production stage, so when β was 10-10 3 , the cumulative oil production at non-Darcy seepage condition was lower than that at Darcy seepage condition, while when β was 10 3 , the results were opposite. Compared with conclusions drawn from the conceptual model, the high-speed non-Darcy showed a more obvious inhibiting effect on water phase, and when the water cut was low, the inhibiting effect of weak high-speed non-Darcy flow on water phase would decrease and oil production might reduce, while when this effect was strengthened, Fig. 6 The effect of high-speed non-Darcy features on production pressure difference oil production would increase and water production would decrease.

Optimization and prediction of injection-production dynamics
S Block is a typical fracture-karst reservoir. The karst caves of different sizes and shapes are developed along the fault zone, and unpacked dense fractures are developed on the two sides (Yun 2021;Zi-Cheng 2020;Bao-Zeng 2020;Fang-Zheng 2018). S Block has been put into trial production since September 2015, and the pressures of some oil wells underwent rapid drop. Hence, it is necessary to optimize the production mode, the energy replenishment timing and   . 8 Single well flow pressure fitting and block cumulative oil production sensitivity analysis the injection/production parameters in block S based on the connectivity prediction model.

(a) Conversion wells selection and well pattern adjustment
Select the connected well groups in block S and set the well spacing to be 0.59-2.21 km to construct the basic injection-production well pattern; Choose the connected units with high conductivity, small controllable reserve and abundant remaining oil. Meanwhile, select the wells experiencing fast flowing pressure drop and with trajectory at the bottom of the reservoir for water injection. As shown in Fig. 9, S-4 and S-22 were selected as conversion wells and three water injection wells INJ_ 1, INJ_ 2, INJ_ 3 were added to improve the basic injection-production well pattern (Fig. 9).
(b) Liquid production Optimization in the injection stage According to Eq. (13) and Table 2, four oil wells in different fault zones were selected for critical production calculation. According to Table 3, the reasonable productivity of S-1 fracture zone was about 50-120t/d.
Take the connected well group as the objective, calculate the reasonable liquid production volume and oil production rate at unit pressure drop, and determine the adjust direction by comparing with the current production performance. As for S-1 fracture-karst reservoir, the production rate should be adjust according to requirements; Based on the above critical liquid production calculation of each single well, and take the liquid production at unit pressure drop as the objective to optimize the liquid production. Results are given in Table 4.
According to Table 4 and Fig. 10, after optimization of production rate, the liquid production of each single well under unit pressure drop increases significantly. Where, the production rates of S-1, S-6 and S-7 increased, S-2, S-5, S-13 and S20 reduced. The liquid production per unit pressure difference in the block S increased by 279.73 m 3 /MPa.
(c) Optimization of water injection parameters 1) Injection timing According to the above derivation of the critical flow pressure at fracture closure, S-1 and S-7 were selected as typical wells. Substitute the rock physical parameters of block S into Eq. (15), and the critical flow pressures at fracture closure in the near well area of S-1 and S-7 were found to be 35.56 MPa and 33.42 MPa, respectively, as shown in Table 5. Therefore, water should be injected when the average formation pressure dropped to 35 MPa.

2) Optimization of injection rate
According to the calculated optimized liquid production rate in Sect. 4.2 the dynamic reserves and the current production status, the interwell remaining recoverable reserves could be calculated. According to Table 6, in this part, a 2-year production prediction simulation was carried out to calculate the oil production at per unit pressure drop (K o ).
Based on the proportion of remaining recoverable reserves and oil production at per unit pressure drop given in Table 7, S-4, S-22, INJ_1, INJ_2, INJ_3 were selected as water injection wells. Moreover, the injection timing when the average formation pressure dropped to 35 MPa was selected and injection/production optimization was conducted based on the prediction of optimized liquid production. A 2-year injection production optimization prediction

(d) Optimization of injection cycle
The main fluid storage spaces in block S are pores and caves of different sizes. The elastic force and the mass force between oil and water in the pores and caves as well as the quality differentiation effect were comprehensively considered. Based on the optimized injection timing and injection production liquid rate, a 64-month cycle water injection prediction optimization design was carried out. In this paper, different cycles including 2 months (1 month of injection and 1 month of well shut in), 4 months (2 months of injection and 2 months of well shut in), 8 months (4 months of injection and 4 months of well shut in), 16 months (8 months of injection and 8 months of well shut in) and 32 months (16 months of injection and 16 months of well shut in) were compared. Determine the water injection cycle by taking the cumulative oil production of the block as the objective. Figure 11 shows that when the injection cycle is set to be 8 months, the cumulative oil production of the block was the maximum. When the cumulative injection volume was the same, the process was a piston type displacement with oil-water interface rising slowly, which was conducive to oil production. When the injection cycle was more than 8 months, the connectivity of the well groups would become better contribute to the overlong continuous water injection and the risks of water channeling and water logging may increase, resulting in low sweep efficiency and low oil recovery. Therefore, the optimal water injection cycle should be 4-8 months.

Conclusions
(1) A connectivity injection/production prediction model for fault-karst reservoirs was established by comprehensively consider the physical characteristics of fault-karst reservoir, the impacts of viscous and inertia forces, and the impacts of high-speed non-Darcy flow. Take the controllable reserves of single well, bottomhole flowing pressure and the cumulative oil production of block S as the fitting indexes, we found that when β was 10 3 -10 4 , the degree of fitting is higher than 85% and the static/dynamic interwell connectivity was consistent with the real case, verifying the accuracy of the proposed model and the connectivity characteristics of connected units were quantitatively characterized.
(2) When β = 0 (Darcy flow), the conductivity value remained unchanged. When β > 0 (high-speed non-Darcy seepage), the conductivity would decrease with time; moreover, non-Darcy seepage became stronger with increasing β value. Meanwhile, for the same medium system (hole-hole, crack-crack), there is no drastic change of flow potential between nodes, and the conductivity decreases linearly. In the different medium system (fracture-cavity), the oil/water interaction produces a drastic change in flow potential, the conductivity decreases nonlinearly, and the change rate decreases with time. In the isolated karst cave, the change rate of conductivity is 0 approximately. (3) When the rate of production liquid was a constant, the high-speed non-Darcy flow might reduce the oil-water flow ability, prolong the water breakthrough time, and significantly affect the water phase seepage ability, therefore, to increase oil recovery and decrease water production. In addition, with the increase in β value, the pressure loss would increase. Thus, the stronger the non-Darcy flow is, the greater the production pressure difference is required. When the bottom-hole flowing pressure was a constant, with the value of β increased, there is a critical value, which changes the law of oil/ water interaction and the energy transfer between injection and production wells. Meanwhile, the high-speed non-Darcy effect mainly occurs in the early and late stage of water drive. (4) The optimal conversion well was selected and three new injection wells were added based on field experience to improve the basic injection-production well pattern and to balance the overall reserve exploitation. In addition, reservoir engineering theory was combined to determine the injection timing, the optimal timing was the time when the formation pressure in block S dropped to 35 MPa. The reasonable productivity of S fracture zone was about 50-120m 3 /d.
(5) Take the constant and stable supply radius as the optimization index and calculate the oil production and water injection at per unit pressure drop the injection rates for S-4, S-22, INJ_1, INJ_2, INJ_3 to be 500 m 3 /d, 350 m 3 /d, 250 m 3 /d, 120 m 3 /d, 300 m 3 /d, respectively. (6) Simulation conducted based on the optimized injection timing. Injection-production rate and dynamic oil saturation change showed that the unit cumulative oil production was the maximum when the injection cycle was 4-8 months.