Oil recovery for fractured-vuggy carbonate reservoirs by cyclic water huff and puff: performance analysis and prediction

Cyclic water huff and puff (CWHP) has proven to be an attractive alternative to improve oil production performance after depletion-drive recovery in fractured-vuggy carbonate reservoirs. However, due to the impact of strong heterogeneity, multiple types of fractured-vuggy medium, poor connectivity, complex flow behaviors and oil-water relationship, CWHP is merely suitable for specific types of natural fractured-vuggy medium, usually causing a great difference in actual oil-yielding effect. It remains a great challenge for accurate evaluation of CWHP adaptability and quantitative prediction of production performance in fractured-vuggy carbonate reservoir, which severely restricts the application of CWHP. For this study, we firstly enable the newly developed fuzzy grey relational analysis to quantify the adaptability of CWHP. With production history of several targeted producers, the accuracy of the proposed method is validated. Based on the traditional percolation theory and waterflood mechanisms in various types of fractured-vuggy medium, a quantitative prediction model for cyclic water cut fwp and increased recovery factor ΔR is presented. The CWHP production performance is discussed by using the Levenberg-Marquardt algorithm for history matching. With a better understanding of the fwp ~ ΔR curve characteristics in different types of fractured-vuggy medium, proper strategies or measures for potential-tapping remaining oil are provided. This methodology can also offer a good basis for engineers and geologists to develop other similar reservoirs with high efficiency.

carbonate reservoirs after water flooding and N 2 flooding by creating specific physical models. Starting from the practice situation of Tahe oilfield, and integrated dynamic evaluation of depletion drive performance, Rong et al. 14 put forward seven major patterns and thirteen sub-classes of remaining oil distribution. The potential-tapping strategies were further presented for various remaining oil distribution patterns. In addition, to tackle the drawbacks of rapid production decline and low recovery rate caused by weak natural energy and sharp water cut rise, many researchers also made great attempts to investigate the feasibility of cyclic water huff and puff (CWHP) 15,16 , N 2 stimulation 17 , N 2 flooding 12,18,19 and N 2 -based foam 20 to enhance oil recovery for fractured-vuggy carbonate reservoirs, as they believed that there should be a certain EOR method suitable for this type of reservoir. The CWHP was proven to be one of the most efficient techniques to improve oil production performance during the late stage of depletion-drive recovery mainly because of its high input-output ratio and relatively strong flexibility. However, pilot tests indicated that, the CWHP was not adapted to all types of fractured-vuggy bodies 10 . With the increase of water injection during CWHP, the oil-yielding effect tends to be worse, and more and more CWHP producers are ineffective, thus causing large amount of remaining oil unexploited at upper positions of producers. Moreover, the coexistence of porous and free-flow domains over a wide range of scales may result in inaccurate description of the complicated fluid flow behaviors in fractured-vuggy media with experimental and theoretical studies in most cases [21][22][23] . How to resolve the adaptability evaluation of the CWHP and understand the actual production performance in various types of fractured-vuggy bodies in order to search suitable potential-tapping strategies or measures remain a great challenge to all of us 24,25 .
In last decades, numerous statistical approaches were proposed to investigate various physical, chemical and biological processes, in which the GRA (Grey Relational Analysis) is treated as an attractive tool to measure the degree of approximation among the data sequences with a dimensionless grey relational grade. The GRA 26 denotes a grey system theory that evaluates the complicated intrinsic relationship between one influential parameter and the others in a specified system. The absolute difference between data sequences is essentially measured to investigate the approximate relationship among sequences with less data. The GRA technique is also employed to understand the synergistic effects of various parameters in order to overcome the shortcomings encountered by other statistical methods 27 . The GRA has a great advantage to determine the intrinsic relationships within less data sequences controlled by various influential factors. Although the GRA technique was successfully used for various areas, e.g., engineering practice [28][29][30] , medicine 31 , hydrocarbon recovery [32][33][34] and microbial production [35][36][37] . To our best knowledge, few studies used this method to study the adaptability of CWHP in fractured-vuggy carbonate reservoirs.
The goal of this study is to understand actual oil recovery in naturally fractured-vuggy carbonate reservoir by combining grey relation analysis of CWHP adaptability and production performance prediction of CWHP

Grey Relational Analysis for CWHP Adaptability
Due to lack of accurate evaluation on the potential of remaining oil after depletion-drive recovery in fractured-vuggy carbonate reservoirs, it is of great difficulty to screen the best potential producers for CWHP. In this work, taking the synergistic effects of drilling data, geological features, dynamic evaluation results and depletion-drive production responses into account, an integrated evaluation index system for CWHP adaptability in fractured-vuggy carbonate reservoirs is firstly developed to improve the decision-making ability of potential tapping. Thereafter, fuzzy grey relational analysis is performed to deal with the problem with high uncertainty.
Evaluation index system of CWHP adaptability. To achieve an accurate evaluation of CWHP adaptability in fractured-vuggy carbonate reservoirs, it is important to consider the synergistic impacts of various static and dynamic data obtained from pilot tests, mainly including the drilling data, geological features, dynamic evaluation results and depletion-drive production responses. A static-dynamic evaluation index system for CWHP adaptability is firstly developed. Table 1 also shows the effect of various evaluation indexes on adaptability of the CWHP. In regard to a targeted fractured-vuggy reservoir with large leakage liquid volume, it is more possible to achieve a higher oil production by cyclic water huff and puff. As the aquifer volume increases, the bottom-water coning during CWHP is usually faster, thus resulting in a more rapid water breakthrough and shorter duration time of primary recovery. Figure 2 shows the curve-fitting procedures for estimating the well-recoverable reserve and aquifer volume to a targeted producer in fractured-vuggy reservoir. Taking the impact of production rate and well bottom-hole pressure into account, we select the Blasingame 38 type curve and flowing material balance (FMB) method 39 for preliminary evaluation. Thereafter, an analytical radial model is developed for single-well production history matching in order to obtain a favorable estimation of well-recoverable reserve and aquifer volume.
To quantify the CWHP adaptability in fractured-vuggy carbonate reservoirs, it is critical to evaluate the impact of various controlling factors and normalize the original data sequences to alleviate the dimension inconsistency. Table 2 displays the ranges of various controlling factors and the normalization of targeted producer's original data sequences. The symbol '↗' denotes to a positive effect of investigated factors on CWHP adaptability, whereas the symbol '↘' shows a negative impact on the CWHP adaptability.
Fuzzy grey relational analysis. Grey relational analysis (GRA) is an efficient statistical technique for determining the degree of approximation among data sequences with a dimensionless grey relational grade. However, Yazdani et al. 40 confirmed that, the usually considered 0.5 of resolving coefficient in the traditional GRA is often inadequate to address the ambiguity found in the available information, as well as the critical fuzziness in decision judgement and preference. It is more suitable to characterize the degree of data certainty using an interval. Therefore, to evaluate the adaptability of CWHP in fractured-vuggy carbonate reservoirs, the newly developed fuzzy grey relational analysis (FGRA) 41 is adopted. The procedures are illustrated as follows.
Step 1: Establish the reference matrix and comparison matrix. The reference matrix is described as where the reference matrix denotes a presumed criterion or optimal values of the investigated factors. If the production performance is better as each factor increases, the maximum value of data sequence is regarded as the reference criterion, otherwise, the reference variable is equal to the minimum value of data sequence. If the dimension of original sequences is m and the CWHP adaptability is investigated under n different influential factors, so the comparison matrix is given by: www.nature.com/scientificreports www.nature.com/scientificreports/ (1) (2)    www.nature.com/scientificreports www.nature.com/scientificreports/ Step 2: Normalize the original data sequences. As the influential factors and the reference values have great dimensions, the following equation should be used for normalization Step 3: Estimate the cosine value of fuzzy membership. For this study, the included-angle cosine method is used, which is not influenced by the linearly proportional data relationship. The similarity of the two factors is achieved according to the value of the included-angle cosine. It is computed as follows: Step 4: Determine the grey relational coefficient ξ ik .
where ψ denotes the resolving coefficient, [0, 1]. The resolving coefficient must satisfy the integrity and anti-interference of the relational coefficient due to that a large or small value is unable to describe the intrinsic relationship of influential factors. The following procedure is used to obtain the favorable resolving coefficient: (a) Compute the mean value of all absolute difference Δ: For this study, if c < 1/3, ψ = 1.25c, else if c ≥ 1/3, ψ = 1.75c.
Step 5: Obtain the Euclidean distance between the reference and comparison matrice, and define the weight of influential factors in the reference matrix using the analytical hierarchy process.
, ] Then, the Euclidean grey relational grade r 2 is obtained by the following formula:  www.nature.com/scientificreports www.nature.com/scientificreports/ Step 6: Determine fuzzy grey relational grade using the fuzzy membership coefficient and the Euclidean grey relational grade, which can be described as Step 7: Using the magnitude of fuzzy grey relational grades, the adaptability of CWHP in different single-well fractured-vuggy units is finally evaluated. Table 3 lists the estimated weight values of influential factors during the evaluation of CWHP adaptability according to analytical hierarchy process (AHP). It indicates that, the dynamic production behaviors during depletion-drive recovery usually have a remarkable impact on the adaptability of CWHP. With respect to a typical fractured-vuggy unit having low water cut before water injection, large well-recoverable reserves and small aquifer volume, it is highly efficient to perform CWHP for further potential tapping.
On the basis of the magnitude of the fuzzy grey relational grades, the adaptability of CWHP in typical single-well fractured-vuggy units are further ranked into four different categories, i.e., excellent, good, poor and bad, denoting to ranges To validate the accuracy of FGRA, the dynamic and static data of three CWHP producers in a fractured-vuggy carbonate reservoir are used to evaluate the adaptability of CWHP. Table 4 lists the values of influential factors and the estimated fuzzy grey relational grades with respect to three targeted producers. The evaluation results are also shown in Fig. 3.
It can be seen from Fig. 3, the larger the area restricted by broken lines, the higher the adaptability of CWHP in single-well fractured-vuggy units. In practice, the cumulative oil production of typical CWHP producers A, B and C are 3.51 × 10 4 t, 3754t and 1433.2t, respectively, which is consistent with the result of fuzzy grey relational analysis, indicating that the static-dynamic evaluation index system and the FGRA method are efficient enough to help decision-making of the best CWHP wells for further potential-tapping of remaining oil in fractured-vuggy carbonate reservoir.   Table 4. Value of influential factors to targeted producers and relevant evaluation results.

Performance Prediction for CWHP
Due to the impact of strong heterogeneity, various types of fractured-vuggy bodies, poor connectivity, multiscale flow behaviors and complicated water-oil relationship, CWHP is merely suitable for specific types of fractured-vuggy bodies, thus causing a great difference in actual production performance. Moreover, with the increasing of water injection during the CWHP, it may suffer from problems of rapid production decline and abrupt increase of water cut, and large amount of remaining oil are still around the top of producers. Therefore, it is essential to quantify the uncertainty and complexity of CWHP production performance in order to explore the proper potential-tapping strategies. For this study, based on the traditional percolation theory and waterflood mechanisms for fractured-vuggy bodies, a quantitative prediction model of cyclic water cut and increased recovery factor is proposed. By using the L-M algorithm for history matching, the actual CWHP production performance for different types of fractured-vuggy medium are further investigated.
Quantitative prediction model for CWHP. With regard to traditional waterflooding reservoirs, the widely used method to describe the relation for oil/water relative permeability ratio and water saturation is the Craft equation, which can be stated as ro rw w where K ro and K rw are the oil-phase and water-phase relative permeability, respectively; S w is the water saturation; a and b are the unknown parameters. However, previous studies indicate that, with respect to ultra-high water cut period, there exists a great deviation between the predicted performance by the traditional Craft empirical model and the observed data 42 , which severely restricts the application of this method. To overcome the inherent drawback, an improved Craft relationship was developed as follows The generalized Darcy's law and water fractional flow equation are described as.
K dp dz a dp dz 1 1 As for single-well fractured-vuggy units, the increased recovery factor by CWHP can be computed as where ΔR is the increased recovery factor; V p is the pore volume, m 3 ; ρ o is the density of crude oil, kg/m 3 ; N is the original oil in place (OOIP), t; B o is the volume coefficient of crude oil; S w0 is the water saturation of fractured-vuggy body during the late stage of primary depletion; S w is the water saturation after CWHP. With regard to a specific fractured-vuggy reservoir, the pore volume, the density of crude oil, the original oil in place (OOIP) and water saturation after depletion-drive recovery are usually constant, and the increased recovery factor is positively proportional to the water saturation after CWHP, i.e., ΔR ∞ S w . Therefore, the quantitative prediction model of cyclic water cut f wp and increased recovery factor ΔR is described in the following.
where a, b, c and d are the unknown model parameters, respectively.
To obtain a satisfactory history matching of the observed production data, a least-square objective function is established, and takes the form of Eq. (19).
where O(m) denotes the objective function; m denotes a (n × 1) vector of the controlling parameters; T is a symbol indicating the transpose of vector or matrix; d obs denotes a (n × 1) vector of observed data; g m ( ) denotes a (n × 1) vector of predicted data; C D denotes the (n × n) covariance matrix of data.
In this paper, optimization is performed using the Levenberg-Marquardt (L-M) algorithm. Besides, a finite difference method is employed to calculate the gradient of objective function at unknown controlling parameters. When using the L-M algorithm to investigate a least-square history matching problems, smooth transmission is usually achieved between the steepest descent algorithm and Newton algorithm. If the least-square objective function is far from the optimal solution, the convergence direction is close to that obtained by the steepest descent algorithm. Otherwise, the convergence direction will be identical to that of the Newton algorithm.
The iteration computation scheme to solve the least-square problem using the L-M algorithm is stated as follows.
where λ is the damping factor; H m ( ) k is the symmetric, half-positive Hessian matrix; I is an identity matrix; δm k+1 is the difference matrix between kth and (k + 1)th iteration.
The following describes how the L-M algorithm is used for optimization of a least-square objection function. First, input the initial damping factor λ 0 . When each iteration is over, it is key to update the damping factor. The basic principle is presented as follows: (1) Calculate the vector of unknown controlling parameters m k+1 . If O(m k+1 ) ≥ O(m k ), the iteration is treated as a failure, and then λ = λ × 10. If O(m k+1 ) < O(m k ), the iteration is successful, and then λ = λ/10. (2) Input the updated damping factor λ into Eq. (20) and perform the next iteration. The above-mentioned iteration is repeated until the termination condition is ultimately satisfied. Figure 4 displays the detailed calculation procedure of analyzing the actual CWHP production performance in fractured-vuggy carbonate reservoirs based on the proposed f wp ~ ΔR quantitative prediction model.

Results and Discussion
Based on production history of CWHP in typical fractured-vuggy carbonate reservoir founded in Tarim basin, the proposed prediction model for f wp ~ ΔR is used to understand the characteristics of cyclic water huff and puff in various types of fractured-vuggy reservoir bodies. Figure 5 shows the ultimate curve fitting effect of typical karst cave. It demonstrates that, the f wp ~ ΔR curve of karst caves seems like a concave profile, denoting to a satisfactory outperformance. Proper strategies are preferable to enlarge the water-free production period and inhibit the sharp bottom-water conning as long as possible. Once water breakthrough is achieved, oil production becomes worse. When the CWHP is ineffective, the nitrogen-based huff-n-puff has proven to be an attractive alternative to displace the attic remaining oil around the top of karst caves.
With respect to eroded pores and caves, the f wp ~ ΔR curve during cyclic water huff and puff obeys a S-shaped or stair-rising law, as shown in Fig. 6, which indicates a relatively great oil production and long water-free production period. Majority of remaining oil unexploited around the top of producers will be recovered after water breakthrough due to the strong bottom water support. During the later stage of production, shut-in to inhibit water coning or water drainage-oil production technique is preferable to improve the CWHP development effect with the greatest effort.
As for the dissolved pores and fissures, the f wp ~ ΔR curve usually conforms to a down-convex law, as shown in Fig. 7. Rapid production decline and low recovery rate are obtained due to weak natural energy and rapid water cut rise. Once water breakthrough is achieved, oil production rate will decrease sharply. The producers for CWHP are under long-time interval production, and the effect of shut-in to control the speed of water conning is negligible in most cases. The unconventional potential-tapping measures, such as large-scale acid fracturing-hydraulic dilation, sidetracking, and nitrogen-based huff-n-puff are recommended for this type of fractured-vuggy medium.  www.nature.com/scientificreports www.nature.com/scientificreports/ conclusions In order to tackle the drawback of inaccurate evaluation on potential of the remaining oil after depletion-drive recovery in fractured-vuggy carbonate reservoirs, we enable the newly developed fuzzy grey relational analysis to quantify the adaptability of CWHP in different types of single-well fractured-vuggy units. Using the production history of several targeted producers, the accuracy of the proposed method is validated, which is of great importance to improve the decision-making ability of potential-tapping.