Multi-scale fracture prediction and characterization method of a fractured carbonate reservoir

The WZ oilfield is characterized by a small production range, low recovery degree, strong reservoir heterogeneity, and complex fracture distribution. At present, there is no quantitative evaluation method for fractures of different scales. This causes problems that include an unclear understanding of reservoirs’ physical properties and remaining oil distribution and seepage characteristics. In this paper, multi-scale fracture prediction and a quantitative characterization method of a fractured carbonate reservoir are studied using three-dimensional seismic imaging logging and regional tectonic stress field distribution data. On the basis of analogs, variance cube, curvature, and the Pompano through-flowline system, large-scale crack recognition is carried out. Combined with the maximum positive curvature attribute, fracture density, and fracture direction interpreted by imaging logging, a small- and medium-sized fracture model is established. Finally, the multi-scale fracture prediction is carried out. This study has important theoretical significance for accurately describing and characterizing the multi-scale fracture distribution law and guiding the fine development of oilfields.


Introduction
The WZ reservoir type is a buried hill, fractured vuggy zone gas-cap block bottom water reservoir. This reservoir type is characterized by strong heterogeneity and complex fracture distribution (Libo et al. 2018). Therefore, accurate description and characterization of multi-scale fracture distribution are of significant importance for future development. In recent years, a large number of researchers have put forward a variety of methods and technologies for the prediction of fractured reservoirs. Dalley proposed that dip and azimuth attributes can be used for seismic interpretation (Dalley et al. 2007) and explained the basic algorithm of these two attributes in detail. Later, these two types of seismic attributes were used by seismic interpretation scholars to interpret small faults, with good results. Curvature attribute is another type of attribute body that is sensitive to discontinuities in seismic data, which can be used to interpret faults and predict fracture tendency and parameters. Bahorich developed coherence cube (Bahorich and Farmer 1995). The core algorithm of the coherence cube calculates the similarity of seismic amplitude volume and highlights dissimilar features to detect the boundary properties. The coherence cube has been further optimized by other scholars, e.g., by an assembly-based coherence algorithm proposed by Marfurt et al. (Marfurt et al. 1998;Gersztenkorn and Marfurt 1999). This algorithm is superior to the original three-channel correlation algorithm and allows the use of lower quality seismic data. Candès presented two fast discrete implementation methods based on the second-generation curvelet transform (Candès et al. 2006). The second-generation curvelet is a new multi-scale analysis tool unconnected to Ridgelet theory (Candès 1999). Characteristically, its implementation process is simple and fast. Ren carried out rock mechanics testing and X-ray computed tomography scanning (Ren et al. 2020) to determine key values for the generation and development of multi-scale fractures in carbonate rocks. Combined with the composite fracture criterion, Ren combined the mechanical model with the numerical simulation software ANSYS to quantitatively simulate the three-dimensional (3D) distribution of multi-scale fracture density and pore size. On the basis of the Roberts algorithm (Roberts 2001), Marfurt (2006) realized the calculation of seismic curvature along the stratum. Subsequently, Al-Dossary and Marfurt proposed the curvature calculation of multispectral volume (Al-Dossary and Marfurt 2006), which can be used to identify small irregular geological features. However, their approach is somewhat complicated. These curvature calculation methods, based on horizontal data, can easily be influenced by picking errors in horizontal interpretation. Klein (Klein et al. 2008) used a similar method to the coherent algorithm with inclination control to calculate the curvature properties. Cook (2016) evaluated Mississippi limestone fractures by correlating volume attributes calculated from surface seismic data with image logs running in six horizontal wells. Pollock (Pollock et al. 2018) proposed a fracture prediction method using common data combined with commercial software and evaluated the fractures of Wolfcamp and Spraberry unconventional reservoirs. Ofoegbu (Ofoegbu and Smart 2019) introduced a program for modeling the initiation, deformation, and propagation of discrete fractures in continuous rock analysis and provided numerical simulation. Nosjean (Nosjean et al. 2020) carried out comprehensive fracture research on a tight reservoir of the Upper Ordovician in Algeria. The study proposed that using 3D-focused diffraction energy technology can greatly assist in characterizing fractures. For the attribute of azimuthal reflected wave amplitude variations, Jin and Stovas (2020), as well as Downton and Roure (2015), discussed its connection with fracture direction, crack density, and infillings. Hunt (Hunt et al. 2010) investigated the accuracy of seismic attributes in predicting fracture density changes in the Nordegg Formation in Midwest Alberta from cores and drilling samples. The study states that the curvature attribute of post-stack seismic data and the azimuthal reflected wave amplitude variations of pre-stack seismic data have strong correlations with imaging log data.
Because of the randomness of reservoir spatial distribution, a single fracture prediction method can only be used for a specific type of formation in the development of a fractured carbonate reservoir. Comprehensive research on limestone reservoir characteristics is lacking, as is targeted prediction methods for fractures of different scales. This research paper employs core, conventional logging, formation micro-imager logging, and geophysical data to identify and quantitatively evaluate fractures in multi-scale limestone reservoirs. The quantitative characterization method of fracture development in a fractured limestone reservoir is established. Concurrently, we further predict the fracture development area of a limestone reservoir and analyze the relationship between fractures and hydrocarbon accumulation. Finally, this study provides theoretical guidance for the further development of the WZ reservoir.

Multi-scale fracture prediction and characterization
There are mainly two aspects for the study of fractures in this reservoir, as shown in Fig. 1. The first is the prediction and acquisition of fracture static parameters. The second aspect is the study of fracture flow characteristics and distribution within the entire area. For large-scale fractures, the seismic attributes suitable for the prediction of large-scale fractures in the study area are optimized by extracting and optimizing the 3D seismic attributes and combining it with the development of faults and the development characteristics of fracture formation in imaging logging. For medium and small-scale fractures, the post-stack 3D seismic data cannot achieve the accuracy of identifying the fractures in this scale. Therefore, spatial fracture development strength is

Fractures with large width and extension
3D seismic data were used to interpret and predict the fault system of the reservoir. On the basis of these data, the development characteristics of large-scale fractures in the study area were identified by multi-attribute volume analysis.

Boundary detection attribute extraction
On the basis of noise reduction, five types of boundary detection attributes, i.e., variance volume, maximum positive curvature, dip offset, similarity, and thinned fault likelihood (TFL), were selected to detect the smoothed amplitude data volume.

Seismic attribute of variance cube
The basic idea of variance cube is using the similarity between adjacent seismic signals to describe the heterogeneity of strata and lithology (Guo 2011). This method is very effective for identifying faults and understanding sand body distribution closely related to reservoir characteristics.
The variance volume can also adopt different timesize windows. In this study, the plane parameter selection is 3, i.e., three seismic data pixels. The longitudinal parameters are 30 and 40, respectively. Figure 2 shows different variance cube attributes extracted from different parameters. By comparing the two images, it can be seen that the sensitivity of variance cube to the selection range of the longitudinal time window is weak. If the parameter of longitudinal smoothing is larger, information of the abnormal body will be lower. This will also make the vertical fault information clearer. After considering the information content and fracture morphology, the longitudinal smoothing parameter was chosen as 30.
To eliminate the influence of formation on the extraction of differential volume, it is necessary to filter by angle. Furthermore, the influence of the in-phase axis in the near horizontal direction on the extraction result of the cube difference is further reduced. In this study, it is reasonable to adopt a 1.5-time window for inclination correction. 2. Curvature seismic attributes Lisle (1994) found that the curvature property was consistent with the open fractures observed from outcrop data. Roberts (2001) described in detail the basic theory of curvature properties and provided the calculation formula for deriving the curvature attribute of the bedding plane. His research shows that the curvature attribute is extremely effective for extracting geometric features such as fault and fracture strike. The physical meaning of curvature is the degree of curvature at a point on the curve. It can be defined in the form of the following differential.
Formula 1 indicates that the calculated curvature value is an absolute value; therefore, there is no posi- tive or negative difference in the mathematical definition of curvature. In geology, however, the curvature is positive and negative. As shown in Fig. 3, when a local layer is a horizontal plane, the normal vectors of the layers are parallel to each other, and the curvature is zero. When the stratum is convex (anticline), the normal vector diverges, and the curvature is positive. When the local stratum is depressed (syncline), the normal vector converges, and the curvature is negative.
To fully express the physical meaning of seismic data, we can define 3D curvature as follows. The plane intersects with the formation surface on a line, and the curvature of all points on the line can be obtained using Formula 1. Therefore, the calculation of 3D curvature is essentially converted to two-dimensional curvature. However, only normal curvature is concerned in the geological field. The phase axis in seismic amplitude data represents the bending and strike of strata; as such, the 3D seismic curvature cube attribute can be obtained by calculating its curvature. In this calculation, the higherorder derivative can be reduced by the following steps: (1) the first derivative can be obtained by suitable initial boundary conditions; (2) the second derivative can be obtained by deriving the first derivative. Because the time window is small, the space surface in this range can be fitted using Formula 2.
In Formula 3, apparent obliquity P (inline) and Q (xline) are calculated by 3D volume scanning. After calculating the values of P and Q, the polynomial coefficients of the quadric surface can be calculated using this formula. The mean curvature and Gaussian curvature can be obtained using Formulas 4 and 5: The most positive curvature and most negative curvature can be obtained using Formulas 6 and 7: The maximum curvature and minimum curvature can be obtained using Formulas 8 and 9: In this study, the maximum positive curvature attribute is used as the research object. When calculating the maximum positive curvature, the vertical window size is 12 and the horizontal line/trace radius is 1. The extracted section (Fig. 4) shows that the vertical discontinuities can be well displayed on the section with the maximum positive curvature attribute. These discontinuities are mainly small faults, and there is less horizontal interference information. This indicates that the curvature attribute can be used to predict faults in this area. Compared with the variance cube, the attributes extracted by the maximum positive curvature are more sensitive to the discontinuous features in the amplitude, and the recognition accuracy is high.

Dip offset attribute
The dip migration attribute can, to some degree, reflect fault information. In this study, we attempted to use this attribute to identify fractures. However, following the extraction, it was found that the information on fault and fracture could not be observed. Accordingly, the method was abandoned (Fig. 5).

Boundary enhancement
The variance cube and curvature cube can identify the large-scale fractures and faults in the study area to some degree. To obtain a better result for fault recognition, it is necessary to enhance the properties of boundary. As shown in Fig. 6, after the 3D boundary enhancement calculation for the two types of attribute bodies, the cal- Fig. 4 Profile of the maximum positive curvature attribute culation of different attributes was found to be quite different. Among them, the effect of variance body was the most obvious. After calculation, the spatial position and continuity of the fault were greatly enhanced under the buried hill surface, which could be easily identified by physical observation. The results of maximum curvature were also improved after boundary enhancement, and the fracture characteristics of space corresponded well with the results of variance volume. After a comprehensive analysis of the pre-processing of seismic data (discussed in the previous paper), we derived a conclusion related to attribute extraction and boundary enhancement; i.e., attribute extraction is an effective method for detecting the inner fracture of a buried hill. However, different boundary detection methods have different effects. By analyzing the slices of variance bodies with different time windows, we found that the sensitivity of variance bodies to the selected range of longitudinal time windows was weak. Additionally, the larger the longitudinal smoothing parameter, the clearer the vertical fault information. Therefore, the variance cube can effectively identify large-scale faults in the buried hill. From the curvature volume section, we can see that the vertical discontinuity can be well displayed on the section with a maximum positive curvature attribute. Moreover, there was also less interference information in the horizontal direction. From the extraction results of the maximum positive curvature and variance cube, we can see that the attributes extracted by the maximum positive curvature are more sensitive to the discontinuous features in the amplitude cube, and the recognition accuracy is high. Therefore, the curvature cube can effectively identify more details in seismic data volume, but it is not as clear as using the variance cube for large fault identification. After extracting the dip migration attribute, we found that the information of fault and fracture could not be observed in the attribute. Accordingly, the dip migration attribute is not suitable for the identification of buried hill reservoir fractures.

Ant extraction
In 1991, Colorni (Colorni et al. 1991) and others first proposed the concept of an ant colony algorithm and abstracted the process as a specific mathematical problem. An ant tracking algorithm is a type of ant colony algorithm that can be used to identify and track continuous information in seismic attributes. This continuous information typically has a linear distribution. The results of the extraction of seismic attributes for concentrated boundary detection were as described in the section "Boundary detection attribute extraction." As noted, the extraction results of seismic attributes from concentrated boundary detection transformed the discontinuous features in amplitude volume into relatively continuous fault information. This information is typically planar and can be sufficiently tracked using an ant tracking algorithm to achieve a 3D visualization of the fault interpretation effect.
The analysis of the geological characteristics of the study area showed that the attributes suitable for fault identification of buried hill reservoirs were the variance cube and maximum positive curvature attributes. The variance volume attribute can identify large-scale fractures more effectively, whereas the maximum positive curvature attribute can identify more specific fracture details. Therefore, in the extraction process, we selected the more sensitive parameter combination in a bid to derive more information. The detailed parameters are shown in Table 1.
In the process of ant extraction, to obtain a clear fracture result, we attempted to use different azimuth filtering ant tracking to originate the initial ant system. As shown in Fig. 7, the strike of faults without azimuth filtering was mainly in an east-west and north-south direction.
To obtain clearer ant system extraction results in these two directions, we filtered the azimuth angles of east-west and north-south directions and extracted the ant body from the square difference body and the maximum positive curvature body (Figs. 8 and 9). In this way, we derived ants with different attributes under different azimuth constraints; we then fused the ants extracted by filtering from different azimuth angles to obtain new ant systems.
After comparing the extraction results of the variancebased and curvature-based ant systems, we found that the variance ant system mainly recognized the large-scale fracture and crack information. These cracks primarily showed two strike groups, near east-west and northeast-southwest. The ant system based on curvature body was able to identify large-scale cracks in the east-west and north-south directions. In the middle area of Fig. 9, the correspondence in the two directions can be observed. This is in good agreement with the fracture formation series on imaging logging.

TFL properties
TFL is a new seismic attribute reflecting the exact location of faults that can generate accurate and sharp sections. Figure 10 shows the use of TFL attributes to extract large-scale fractures. According to the results of the TFL attribute extraction, the attribute can identify two groups of large-scale fractures in the near-east-west and nearnortheast-southwest directions. In the southern part of the study area, due to the influence of the primary control faults, the large-scale fractures were mainly east-west, accompanied by some small-scale north-south fractures. In the northern part of the study area, there were mainly some large-scale fractures in the north-south direction. On the basis of the rose diagram of fractures interpreted by A1h and A2h imaging logging, it can be seen that the results of TFL attributes show good correspondence with imaging logging.
From the spatial distribution of large-scale fractures, it can be seen that two groups of fractures developed near the A1h and A2h horizontal wells, further indicating that the  . 7 The results of ant tracking extraction based on variance volume attributes without azimuth filtering Fig. 8 The final result of the ant system extraction based on variance body Fig. 9 The final result of the ant system extraction based on curvature body fracture prediction results of the TFL attribute showed good agreement with imaging logging (Fig. 11). After comparing the boundary attribute extraction, the ant system, and the TFL attributes to provide sharpened fault maps on horizontal or vertical slices, we found that TFL produced seismic attributes with sharp edges, which is extremely suitable for fault system interpretation. Accordingly, TFL attributes represent an effective method for predicting large-scale faults in the study area. Finally, TFL attributes can be selected to predict and study large-scale fractures.

Calculation of curvature relative strength of cracks in space
OpendTect software was used to transform the curvature attribute into relative crack density. The corresponding situation for curvature volume generated by seismic data and faults in the study area is shown in Fig. 12. The spatial distribution of crack strength and the azimuth angle of the curvature body are shown in Fig. 13. The obtained fracture development zones were essentially concentrated near faults and were consistent with the trend of large fault systems, which fully conformed to the development law of carbonate fractures (Ren et al. 2020).

Single well fracture description and fracture density correction
The imaging logging data of two horizontal wells in the study area are used to describe the fracture development on a single well (Fig. 14 for logging data and Fig. 15 for fracture development). First, we created statistics on the development position, azimuth, and dip angle of each fracture in a single well. We then found that there were two groups of fractures in the A1h well, in a near-north-south direction and a near-east-west direction. The fracture azimuth angle of the A2h well was similar to that of the A1h well, but the fracture azimuth angle of the A2h well was approximately 20° anticlockwise from that of the A1h well.
Generally speaking, two groups of fractures primarily developed. In the first group, the azimuth of fractures is near-north-east to south-south-west, which is perpendicular to the direction of the main fault and retains a good corresponding relationship with large-scale fractures. In the second group, the fractures are nearly east-west trending, where the A2h well has anticlockwise deflection. This group of fractures is along the fault direction; as such, it may be the fracture associated with the fault.

Calculation of crack density in three-dimensional space
Through the development of fractures on imaging logging data, we can calculate the development intensity of fractures in a single well. The calculated strength is the number of fractures per meter, defined as P 10 . In the section "Calculation of curvature relative strength of cracks in space,"  the relative strength of fracture development was obtained on the basis of the curvature attribute. By comparing the fracture strength of a single well with the fracture density calculated by imaging logging, we can establish the relationship between them and then calculate the fracture density in space using this relationship. Figure 16 shows that the relationship between the two is linear. This relationship can be given as follows: By applying relationship 10 to 3D space, we can obtain the fracture density in 3D space; however, this fracture density represents linear density. When we model fractures, bulk density (P 32 ) is typically needed. Therefore, Eq. (11) is used to convert the fracture linear density into the fracture volume (10) P 10 = 9.15 × P curv + 0.21 density, that is, the total fracture area per unit volume. Furthermore, C 31 is the conversion coefficient from linear density to bulk density.
In the section "TFL properties," we found that the largescale fracture strike extracted by the TFL attribute had a good corresponding relationship with the fracture trend on imaging logging. Therefore, the TFL attribute was used to fit the azimuth angle of two groups of fractures in space, which was subsequently used for fracture modeling (Figs. 17, 18  and 19).
On the basis of the discrete fracture modeling function in Petrel software, combined with the previous fracture description and statistical parameters, we established two groups of fractures, near east-west, and near north-south.
(11) P 32 = C 31 P 10   Tables 2 and 3. The final small-and medium-sized discrete fracture model in the study area is shown in Fig. 20.

Conclusion
1. For fractured reservoirs in a buried hill, attribute extraction is an effective method to detect internal faults. Different boundary detection methods have diverse effects; variance can effectively identify large-scale faults in a buried hill; curvature cube can effectively identify more  Correspondence between the fracture density of the A1h well, the A2h imaging log, and the final fracture body density The dip migration attribute is not suitable for fracture identification of a buried hill reservoir. 2. For the reservoir in the study area, the ant system of the variance cube mainly identified the large-scale fracture and fracture information. These fractures primarily showed two groups of strike, near east-west, and northeast-southwest. The ant system based on curvature cube was able to provide good identification of the large-scale cracks in the east-west and north-south directions. 3. After comparing and optimizing many large-scale fracture identification methods, such as square difference ant body, curvature ant body, and TFL, we found that TFL attributes provided sharpened fault maps on horizontal and vertical slices. Furthermore, TFL produced seismic attribute bodies with sharp edges, which is well suited for fault system interpretation. Accordingly, it is reliable to identify large-scale fractures using the TFL attribute volume in this study area. 4. Using the post-stack inversion data to calculate the maximum positive curvature attribute, we effectively identified the small-and medium-sized fractures by fitting the strength and azimuth attributes.
Funding The work in this paper is supported by the Grant from the National Natural Science Foundation of China (No. U1762107) and Science and Technology Program of Sichuan Province, China (No. 2019YJ0425).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.