Before we go to the details of this study, a brief review of the steering pyramid method is given here. The steering pyramid method is an algorithm for multi-scale, multi-directional decomposition and reconstruction of seismic data. The method decomposes an image into a group of images (called Levels) of the same size as the original one but with different resolution or different frequency band (Mathewson and Hale, 2008). In the process of decomposition, different radial filters and directional filters are designed to enhance data intrinsic features of discontinuity as well as the original structure of the input image remain unchanged. The following filter formulas are defined based on different angles: G k = ∑ m ∑ n W k G k − 1 2 x + m , 2 y + n (1) f θ x ,   y = ∑ j = 1 M k j θ f θ j x ,   y (2) where G k represents the k _th level of the pyramid, x , y represents an array of image pixels from different levels, and W k represents a window function with low-pass characteristics. f θ x ,   y is a function of the controllable filter in the θ direction, which can be obtained by combining the interpolation function k j θ in the θ direction and the linear f θ j x ,   y of the basic function in the θ direction, and M is the number of the basic functions required for steering. k j θ are the weighting functions.

In the seismic application of the decomposition method, seismic data are normally decomposed into 7 levels of different frequency band. The levels holding high-frequency can be stacked to one, and for example Level 0 and 1 are stacked to form a high-frequency stack, which improves seismic data resolution and enhances the discontinuity for the purpose of identification of edges, joints, faults, fractures or cavities. Zhao et al. (2021) used the steerable pyramid technology on seismic data to recognize geological structure. This method effectively removes background noise and enhances the potential geological structures.

Following the seismic decomposition and level stack, a comparison analysis is performed to evaluate the benefits of the process in both frequency domain and time domain, depending on the study objective. For example, high-resolution stacks can exhibit clearer or more consistent images of fault development in 3D space for fault interpretation. This process is shown in a workflow of Figure 1.

Once the seismic decomposition is complete, seismic attributes can be extracted from the stack with high resolution. For fracture identification in this paper, the seismic attributes, variance, curvature and mean consistent curvature are selected.

The variance attribute is a seismic attribute that describes the variability of seismic amplitudes in neighboring seismic traces. It is more sensitive to fractures caused by faults with strong discontinuities and is suitable for identifying faults at median-scale (Zhang et al., 2021). Curvature can provide information about the spatial distribution of structural features in a rock formation, and is most sensitive to fractures developed by fold; it is suitable for identifying large-scale fractures (Al-Dossary and Marfurt, 2006). The mean consistent curvature, a specific curvature derivative, is a type of seismic attribute that measures the curvature of a seismic horizon or event in a consistent manner throughout the entire seismic volume. The mean consistent curvature attribute is designed to be more sensitive to small-scale fractures and other subtle structural features, as it is able to capture variations in the seismic response that are not evident in the original seismic data (Chilès and Marsily, 1993).

The attributes above represents different characteristics of fractures or fracture networks (e.g., geometry, scales and accuracy). To combine the advantages of each of the seismic attributes, a composite seismic attribute with high accuracy is generated using a linear superposition method. The steps to achieve this are as follows.

① Seismic attributes are normalized to ensure that they be superposed.

Assuming that the geological survey large fault area S is 20 km^2. The curvature attribute calculation obtains the fault area S 1 as 10km², and the variance attribute calculation obtains the area S 2 as 18 km². Then the correlation r 1 between the curvature attribute and the fault is 0.5, and the correlation r 2 between the variance and the fault is 0.9. From the above equations, it can be deduced that the weight ω 1 of the curvature attribute is 0.357, and that of the variance ω 2 is 0.643.

In this section, the composite seismic attribute generated from the previous section is utilized as the fracture location constraint in the DFN modelling (Chopra and Marfurt, 2013) by introducing the Poisson process. DFN modelling is a numerical modelling approach used to simulate fluid flow and transport in fractured rock masses.

The principle of the Poisson process is a mathematical model used to describe the occurrence of random events over time or space. It is particularly useful in modelling systems where events occur randomly, or in this paper, fractures distribute randomly. Therefore, the Poisson process can be used in DFN modelling to generate random fractures in a rock mass to create a DFN model simulate the development of natural fractures or fracture network in rocks.

When generating fractures in 3D space, the flow chart of the traditional approach is shown in Figure 2. In this paper, we use the composite seismic attribute to constraint fracture location distribution as the attribute is much related to the development of fractures or fracture network. The steps are as follows.

① The composite seismic attribute is regularized, shown in Eq. 4. The generation probability P r of fractures per unit volume can be calculated as,

Overall, using the Poisson process as an input for DFN simulations can help create more realistic and representative models of fractured rock masses, which can be useful for understanding and predicting fluid flow and transport in geological systems.

The study area, the Bonan Depression, is located in the middle of the Zhanhua Oil field, in East China. In spite of the long production history, the geological structure and fault systems are not fully understood due to the presence of the complex faults and fault system. Seismic interpretation shows that the depression is crossed by two major faults indicated by two bold red curves (Kang et al., 2002; Wang and Zhang, 2023) shown in Figure 3. Many minor faults develop as branches of the two majors.

The study area is approximately 16.2 km from east to west and 9.2 km from north to south, and the thickness of the target layer, the Shasi Interval, is around 15 m. High-angle tectonic fractures are developed within the distance of hundreds of meters of the faults, and there are six wells in the study area. The reservoir model is not effective based on the traditional DFN modelling approach, and the workflow proposed in the previous section will be tested and verified in the study area in order to improve the understanding of the reservoir model.

For the case application, seismic data from the Bonan Depression area were decomposed, and seismic attributes were extracted before the composite attribute was generated. The steerable pyramid technique decomposed the seismic data into seven levels, from Level 0 to 6. Level 0 contained much noise and was not considered. The higher levels (Level 4, 5, and 6), holding low-frequency content, indicated the large-scale trend of the fault system. The other levels (Level 1, 2, and 3) of the high-frequency band were stacked to improve the resolution and prediction of small-scale faults. To examine the workflow, crossline 803 was extracted and presented in Figure 4A. The section stacked with Level 1, 2, and 3 is shown in Figure 4B for comparison. The area highlighted in the ellipses indicates the improved clarity of the two faults (enhancement of the discontinuity), which facilitates seismic interpretation.

Following the decomposition, three seismic attributes, variance, curvature and mean consistent curvature are extracted using the stack. Figure 5B shows the variance for the window parameter of 3*3. Figure 5D shows the variancefor a vertical radius parameter of 12. Figure 5F shows the mean consistent curvature for the window parameter of 3*3. Figures 5A, C, E are the three attributes extracted from the original seismic using the same parameters. It can be seen that the attributes from the original seismic are much noisier, compared with the ones from the stack, from which the outlines of the faults and fault system are more clearly observed. Since fractures develop along the faults or within the fault system. These clear images will reflect more logical fracture distribution if they are utilized to constrain the DFN modelling in this area.

With the extracted seismic attribute, a composite attribute is computed according to the Eq. 3 in the previous section shown in Figure 5G. The composite attribute combines the advantage of the three, reflecting the potential fractures along the fault or within the fault system. And it will be used as the final resulting attribute constrain the DFN modelling in the next section.

With the composite attribute extracted, it is served as the location constraint of the DFN modelling in this section. The modelling process is constructed with the approach presented in the section II(C). In addition, Fracture density distribution (or the number of fractures in one location) is created from wells in the method presented in (Jian et al., 2021) and the details are not discussed here as we focus on the location constraint. By examining the fracture properties at the wells, its orientation is determined as 30° and 135° with the variance of 0.1. The fracture dip angle is mainly high-angle, and the mean value and variance are 75° and 0.5, respectively. In the process of fracture generating, each fracture is created with 4 edges, and the mean fracture length is 10 m, and its width is around 3–4 m.

With the parameters above, the DFN model created with the location constraint is shown in Figure 6B. Figure 6A is the model built without the constraint. By comparison, fractures in Figure 6A distribute across the full area; however, in tectonic sense, they should only emerge along the faults or within the fault system. Therefore, with the location constraint, Figure 6B indicates the consistency of the spatial relation between fractures and the regional faults. For example, many fractures develop in the two ellipses (see Figure 6A) there are no faults, which is not correct in geological sense; at the same locations (the two ellipses in Figure 6B), few fractures develop as the location constraint is applied.

As the same time, we compare the number of fractures near the well with the daily fluid production. The number of fractures within a square with the side length of 100 m, centered at the well. The number of fractures for Wells Y1–Y6 are 2, 4, 21, 14, 14 and 9, respectively. The corresponding daily fluid productions per well were 1.18, 5.19, 27.9, 28.77, 13.2 and 18.5. The fracture distribution at the wells is consistent with the production with the correction coefficient of 0.755 (Figure 6C). This demonstrates the effectiveness of the DFN method.