Discrete Fracture Network Modelling in a Naturally Fractured Carbonate Reservoir in the Jingbei Oilfield , China

This paper presents an integrated approach of discrete fracture network modelling for a naturally fractured buried-hill carbonate reservoir in the Jingbei Oilfield by using a 3D seismic survey, conventional well logs, and core data. The ant tracking attribute, extracted from 3D seismic data, is used to detect the faults and large-scale fractures. Fracture density and dip angle are evaluated by observing drilling cores of seven wells. The fracture density distribution in spatiality was predicted in four steps; firstly, the ant tracking attribute was extracted as a geophysical log; then an artificial neural network model was built by relating the fracture density with logs, e.g., acoustic, gamma ray, compensated neutron, density, and ant tracking; then 3D distribution models of acoustic, gamma ray, compensated neutron and density were generated by using a Gaussian random function simulation; and, finally, the fracture density distribution in 3D was predicted by using the generated artificial neural network model. Then, different methods were used to build the discrete fracture network model for different types of fractures of which large-scale fractures were modelled deterministically and small-scale fractures were modelled stochastically. The results show that the workflow presented in this study is effective for building discrete fracture network models for naturally fractured reservoirs.


Introduction
Fracture prediction is important for evaluating and developing fractured reservoirs [1].However, it is difficult to characterise fractures properly because of their extremely heterogeneity and complexity features [2,3].Geologists and engineers tried to use different types of data to describe the fractures, e.g., seismic attributes, geophysical logs, drilling cores, and outcrops [4,5].A seismic edge detection method was used to conduct stratigraphic and structural interpretations of geological features [6,7].However, some subtle faults and fractures may cause minor changes in the seismic waveform which are not easily correlated using the conventional interpretation of the seismic section.Thus, many seismic interpreters tried to use seismic attributes.Dip and azimuth attributes were put forward by Dalley et al. [8], who described the basic algorithm for using dip and azimuth attributes to interpret the faults.Later, dip and azimuth attributes were used in identifying subtle faults, as well [6,9].More recently, curvature attribute was widely used in delineating faults [6,10,11] and predicting fracture orientations and distributions [7,12,13].Coherency attribute developed by Bahorich and Farmer [14], can calculate the similarity of the seismic amplitude volume to enhance the seismic edge.Marfurt et al. [15] and Gresztenkorn and Marfurt [16] developed an algorithm called semblance-based coherency, which allowed us to analyse data of lesser quality than our original three-trace cross-correlation-based algorithm.An ant tracking attribute is based on the ant colony algorithm which mimics the behaviour of ants to track the continuous feature [6].This algorithm is powerful when it is used to interpret subtle faults and fractures [4,6,[17][18][19][20].In recent years, some researchers tried to combine multiple attributes, which were used to detect the seismic edges, to enhance the interpretation for the subtle faults and fractures.Nasseri et al. [4] used chaos, variance, dip deviation, and coherency to detect the discontinuities of the amplitude volume which then were used to interpret the subtle faults and fractures; then the ant colony algorithm and fuzzy c-means clustering were used to enhance these discontinuous features (subtle faults and fractures).Basir et al. [6] employed dip, similarity, and curvature attributes to detect the discontinuity of the seismic volume and then used the ant tracking approach to enhance the seismic edge.Final faults were interpreted by combining these attributes through artificial neural network (ANN) algorithm [21][22][23][24].
Seismic attributes cannot identify small-scale fractures (SSF) because of the resolution restriction of seismic data volume [25].However, borehole data, e.g., wireline logs, drilling and core data, provide information of fractures directly or indirectly.Therefore, researchers studied the SSFs by using core observation, image log [26][27][28] and geophysical logging interpretation [29][30][31] combining some mathematical methods [32].Image log interpretation is an advanced technique to detect and evaluate the fracture intensity and orientation.However, there is no image log in this study area, which forces us to assess the fractures using conventional geophysical logs.Geophysical logs, e.g., temperature log, caliper (CAL), spontaneous potential (SP) curve, natural radioactivity logs, density logs (DEN), compensated neutron logs (CN), and acoustic logs (AC), etc., can reflect fracture features [33].
Fracture density is an essential parameter to characterise the fracture occurrence.Martinez [34] provided a fuzzy logic technique to predict the fracture density using conventional geophysical logs.Al-Anazi and Babadagli [29] predicted the fracture density using ANN model of which static geophysical logs and dynamic production data were combined together.Tokhmchi et al. [30] predicted the fracture density by converting the conventional geophysical logs, e.g., CAL, AC, DEN, and photoelectric (PHE), to energy logs; then a linear relationship between energy logs and fracture density was established and then used to predict the fracture density.Zazoun [31] predicted the fracture density for un-cored wells by using wireline logs, e.g., AC, gamma ray (GR), DEN, CN, CAL and depth logs, and measured fracture density data from core observations.All of these methods predict the fracture density successfully and quantitatively at the depth point.Nevertheless, it is still a challenge and shortage in the literature to predict the fracture density distribution in 3D.Geostatistics is a widely used method in analysing the trend and structure of the spatially distributed data.From the controlling factors of this trend and structure, a geologically reasonable interpolation scheme can be used to build 2D/3D geological models [35,36].In this study, we tried using geostatistics to predict the spatial distribution of the fracture density by combining ANN algorithm.
Fracture modelling is to simulate the fractures and the fluid flow in the fracture systems [5].Two models, e.g., equivalent continuous model and the discrete fracture network (DFN) model, were used to describe the fracture distribution.The equivalent continuous model was brought forward by Barenblatt and Zheltov [37] and developed by later researchers [38][39][40].In this method, most fluid stores in the matrix and flow through fractures; the fluid flow between matrix and fractures is presented by a transfer function [41].DFN modelling is an advanced approach for fracture modelling, of which different sets of fractures could be created.Each fracture is presented by a plane with specific parameters such as dip angle, dip azimuth and aperture, etc.The fracture model built by DFN approximates the real fracture distribution in plays [42,43].Both of the deterministic and stochastic fractures [44,45] can be generated using DFN.The petrophysical properties, e.g., porosity and permeability, can be calculated by upscaling the DFN model [27,46,47] into cells in order to simulate the reservoir fluid flow [48].In this study, we tried to use seismic data to detect and extract the subtle faults and large-scale fractures because the fluid flow in the study area is mainly through faults and fractures [49,50].Geophysical logs and core observations were used to predict the parameters of small-scale fractures.The fracture density is calculated by an ANN model based on the spatially distributed geophysical properties.Figure 1 shows the workflow.properties, e.g., porosity and permeability, can be calculated by upscaling the DFN model [27,46,47] into cells in order to simulate the reservoir fluid flow [48].
In this study, we tried to use seismic data to detect and extract the subtle faults and large-scale fractures because the fluid flow in the study area is mainly through faults and fractures [49,50].Geophysical logs and core observations were used to predict the parameters of small-scale fractures.The fracture density is calculated by an ANN model based on the spatially distributed geophysical properties.Figure 1 shows the workflow.

Geological Settings for Jingbei Oilfield
The Jingbei Oilfield, with a proved area of 14.1 km 2 for oil, is located in the north of the Jinganbao tectonic belt in the Damingtun Depression, Northeast China (Figure 2a).The buried depth of the oil zone, which belongs to the Proterozoic formation [51] and partly to the Archaeozoic formation, ranges from 2450 m to 3100 m [49,52].Until 2011, 7.57 Mt of oil had been produced.The lithologies for the Proterozoic formation are quartzite, dolomite, slate, and limestone while the Archaeozoic formation beneath consists of mixed granite and metamorphic rocks (Figure 2b).Fractures and faults are widely distributed in the reservoir because of the multiple tectonic events [53].Fractures are the main space for oil in the reservoir and tectonic fractures are the dominant fracture type [52].

Geological Settings for Jingbei Oilfield
The Jingbei Oilfield, with a proved area of 14.1 km 2 for oil, is located in the north of the Jinganbao tectonic belt in the Damingtun Depression, Northeast China (Figure 2a).The buried depth of the oil zone, which belongs to the Proterozoic formation [51] and partly to the Archaeozoic formation, ranges from 2450 m to 3100 m [49,52].Until 2011, 7.57 Mt of oil had been produced.The lithologies for the Proterozoic formation are quartzite, dolomite, slate, and limestone while the Archaeozoic formation beneath consists of mixed granite and metamorphic rocks (Figure 2b).Fractures and faults are widely distributed in the reservoir because of the multiple tectonic events [53].Fractures are the main space for oil in the reservoir and tectonic fractures are the dominant fracture type [52].
The Proterozoic formation undertook multiple tectonic movements [50,53].At the early stage of the Paleozoic, main tectonic movements were upright, which lead to the erosion of the Palaeozoic formation.During the Indo-Chinese epoch, some wide N-S folds were formed in the Proterozoic formation, due to the bearing of N-S directional pressure.In Yanshanian, the study area suffered NW-SE directional compressive stress because of the collision of the Kula-Pacific plate compressing toward the west that led to some new N-E folds.In the early Himalaya epoch, a series of main faults with a near N-E direction were generated under the effect of tensile stress which caused by the Pacific plate diving beneath the Eurasian plate.The Proterozoic formation undertook multiple tectonic movements [50,53].At the early stage of the Paleozoic, main tectonic movements were upright, which lead to the erosion of the Palaeozoic formation.During the Indo-Chinese epoch, some wide N-S folds were formed in the Proterozoic formation, due to the bearing of N-S directional pressure.In Yanshanian, the study area suffered NW-SE directional compressive stress because of the collision of the Kula-Pacific plate compressing toward the west that led to some new N-E folds.In the early Himalaya epoch, a series of main faults with a near N-E direction were generated under the effect of tensile stress which caused by the Pacific plate diving beneath the Eurasian plate.

Data
A 3D seismic survey with an area of 60 km 2 was collected from the oil company.The average dominant frequency of data was 30 Hz.The seismic data were used to detect subtle faults and large-scale fractures.More than 130 wells are drilled in the studied area.Geophysical logs of AC, DEN, CN, and GR from 110 wells were collected and used in this study.The 223 m cores from seven cored wells (locations are shown in Figure 2a) were observed and analysed.Geophysical logs and core observation data were used to identify the features, e.g., fracture density and fracture dip angle of small-scale fractures at the borehole scale.

Condition of the Seismic Data Volume
Seismic attributes are often sensitive to noise, especially for the buried-hill reservoir of which the seismic data quality is severely affected by the erosion surface.Therefore, running a spatial filtering is imperative to reduce the noise in the signal.In this study, we used a structural smoothing approach to condition the original seismic volume.The structural smoothing method smooths the input signal guided by the local structure to increase the continuity of the seismic reflectors [54].Table 1 lists all parameters in the seismic volume condition.Principal component dip and azimuth computation were used to determine the local structure.Gaussian smoothing was then applied

Data
A 3D seismic survey with an area of 60 km 2 was collected from the oil company.The average dominant frequency of data was 30 Hz.The seismic data were used to detect subtle faults and large-scale fractures.More than 130 wells are drilled in the studied area.Geophysical logs of AC, DEN, CN, and GR from 110 wells were collected and used in this study.The 223 m cores from seven cored wells (locations are shown in Figure 2a) were observed and analysed.Geophysical logs and core observation data were used to identify the features, e.g., fracture density and fracture dip angle of small-scale fractures at the borehole scale.

Condition of the Seismic Data Volume
Seismic attributes are often sensitive to noise, especially for the buried-hill reservoir of which the seismic data quality is severely affected by the erosion surface.Therefore, running a spatial filtering is imperative to reduce the noise in the signal.In this study, we used a structural smoothing approach to condition the original seismic volume.The structural smoothing method smooths the input signal guided by the local structure to increase the continuity of the seismic reflectors [54].Table 1 lists all parameters in the seismic volume condition.Principal component dip and azimuth computation were used to determine the local structure.Gaussian smoothing was then applied parallel to the orientation of this structure [55].Then, we enhanced the seismic spatial discontinuities by calculating seismic variance which is an edge detecting method.The edge means the lateral discontinuities of seismic amplitude.Ant tracking [6] is used to detect and extract the continuous features from the seismic attributes.The ant-tracking algorithm is based on the ant-colony algorithm to capture trends in noisy data [6].Non-structural features, such as noise and channels are less likely to be captured by the ant-tracking algorithm because these features usually have internally chaotic textures which are not continuous.
Several parameters, e.g., the initial ant boundary, ant track deviation, ant step size, illegal steps allowed, legal steps required, and stop criteria, significantly affect the extraction result by the ant tracking algorithm [55].The initial ant boundary controls how closely the initial ant agents can be placed within the volume.Ant track deviation defines the searching tolerance of ant on their side of its tracking direction.A larger ant track deviation value allows more connections for the ant agent.The step size, normally using 3, is how far an ant advances for each increment of its search.Illegal steps allow the ant to search beyond its current location when an edge has not been detected.The legal steps required parameter works in combination with the illegal step parameter, by requiring the selected number of valid steps after an illegal step.The stop criteria parameter also works in combination with the illegal steps, and is used to terminate an ant's advance when too many illegal steps have been taken.
In the software, the parameters' combination of the ant tracking algorithm has three options, known as passive, custom, and aggressive parameters [55].As the studied reservoir is a buried-hill reservoir, the reflection is weakened by the erosion surface.Therefore, the aggressive parameters shown in Table 2 were selected for detecting the discontinuity in ant tracking processes.The aggressive option allows the ant agents more freedom to explore and detects more subtle connections.

Extraction of the Faults and Fractures
Automatic fault-extraction technology can reduce human intervention and improve accuracy and efficiency in fault interpretation [56].However, the result of automatic extraction is restricted by the parameters, such as extraction sampling distance, extraction sampling threshold, extraction background threshold, deviation from a plane, connectivity constraint, minimum patch size, and patch down sampling.
Extraction sampling distance sets the minimum distance between extraction seed points.Distance is measured in voxels and defines a radius around the extraction point.The extraction sampling threshold sets the minimum signal level from which to create extraction points.The values range in percentage for signal minimum and maximum.The extraction background threshold sets the minimum signal level to be incorporated into fault estimate.Deviation from a plane controls how far a fault may deviate from a plane surface fit to the data.The value defines the number of voxels to be searched on either side of the orientation estimation plane.For minimum patch size (points), fault patches which contain fewer than this number of points will be excluded from the fault patch set.
Patch down sampling (voxels) controls the density of points within each patch [55].Table 3 lists the parameters which were used in extracting of the fault patches in this study.Artificial neural networks (ANNs) are computer models that mimic the functions of the human nervous system through some parallel structures comprised of non-linear processing nodes which are connected by weights [57].These weights establish a relationship between the input and output of each node in the ANNs [58].These systems process the data and then learn the relationships between the given data in a parallel and distributed pattern.Hence, ANN is robust in capturing the complex relationships among different parameters [31].ANN was also used in seismic data processing and interpretation, geophysical logging interpretation, reservoir mapping, and engineering [58].
In order to control the spatial distribution of fracture density, the ant-tracking attribute log was extracted from the ant tracking attribute as a synthetic wireline log.The ANN model was then built with geophysical logs of AC, DEN, CN, GR, and ant attribute logs as the inputs and fracture density as the output (see Figure 3).Note that the ANN algorithm is an embedded tool in Petrel TM (Schlumberger, Houston, TX, USA) in which the neuron number in the hidden layer is determined from trial and error.The parameters of a maximum number of iterations, error limit, and cross-validation percentage of samples in building the supervised ANN model were set as 20%, 10%, and 50%, respectively.The selection of 50% of samples for cross-validation is used to prevent overtraining/overfitting [59].
Energies 2017, 10, 183 6 of 19 (points), fault patches which contain fewer than this number of points will be excluded from the fault patch set.Patch down sampling (voxels) controls the density of points within each patch [55].Table 3 lists the parameters which were used in extracting of the fault patches in this study.Artificial neural networks (ANNs) are computer models that mimic the functions of the human nervous system through some parallel structures comprised of non-linear processing nodes which are connected by weights [57].These weights establish a relationship between the input and output of each node in the ANNs [58].These systems process the data and then learn the relationships between the given data in a parallel and distributed pattern.Hence, ANN is robust in capturing the complex relationships among different parameters [31].ANN was also used in seismic data processing and interpretation, geophysical logging interpretation, reservoir mapping, and engineering [58].
In order to control the spatial distribution of fracture density, the ant-tracking attribute log was extracted from the ant tracking attribute as a synthetic wireline log.The ANN model was then built with geophysical logs of AC, DEN, CN, GR, and ant attribute logs as the inputs and fracture density as the output (see Figure 3).Note that the ANN algorithm is an embedded tool in Petrel TM (Schlumberger, Houston, TX, USA) in which the neuron number in the hidden layer is determined from trial and error.The parameters of a maximum number of iterations, error limit, and cross-validation percentage of samples in building the supervised ANN model were set as 20%, 10%, and 50%, respectively.The selection of 50% of samples for cross-validation is used to prevent overtraining/overfitting [59].

Discrete Fracture Modelling
In this study, two ways were used to generate the DFN model.For the fractures in which the length is longer than 200 m, a deterministic method was used to generate their distribution by using the extracted discontinuity from the ant tracking attribute, while a stochastic method was adopted to model the fractures shorter than 200 m.The threshold selection of 200 m is explained in Section 3.3.In the stochastic process, the DFN models are crucially affected by the parameters, e.g., fracture density, fracture length, fracture dip angle, and fracture azimuth.

Discrete Fracture Modelling
In this study, two ways were used to generate the DFN model.For the fractures in which the length is longer than 200 m, a deterministic method was used to generate their distribution by using the extracted discontinuity from the ant tracking attribute, while a stochastic method was adopted to model the fractures shorter than 200 m.The threshold selection of 200 m is explained in Section 3.3.In the stochastic process, the DFN models are crucially affected by the parameters, e.g., fracture density, fracture length, fracture dip angle, and fracture azimuth.In the core observation process, we counted all of the fractures occurring on the surface of the cored rock; meanwhile, the dip angle of each fracture is measured.Fracture density on the core surface was then calculated using the area of the fractures of each core divided by the volume of the core which is equal to the core length times the borehole's area.

Seismic Conditions
The original amplitude volume of the buried hill is severely affected by the erosional surface so that the amplitude continuity is interrupted (Figure 4a).Hence, the structure smoothing process was used to reconstruct this continuity by using the Gaussian smoothing algorithm (Figure 4b).Then the variance algorithm was used to detect the edge which represents the discontinuity of the amplitude volume (Figure 4c).In the time slice of the variance attribute volume, we can see there are many discontinuities in the buried hill zone (black zones in Figure 4c).Finally, the 3D edge enhancement method was used to enhance the edge detected in the variance attribute (Figure 4d).

Core Observation and Analysis
In the core observation process, we counted all of the fractures occurring on the surface of the cored rock; meanwhile, the dip angle of each fracture is measured.Fracture density on the core surface was then calculated using the area of the fractures of each core divided by the volume of the core which is equal to the core length times the borehole's area.

Seismic Conditions
The original amplitude volume of the buried hill is severely affected by the erosional surface so that the amplitude continuity is interrupted (Figure 4a).Hence, the structure smoothing process was used to reconstruct this continuity by using the Gaussian smoothing algorithm (Figure 4b).Then the variance algorithm was used to detect the edge which represents the discontinuity of the amplitude volume (Figure 4c).In the time slice of the variance attribute volume, we can see there are many discontinuities in the buried hill zone (black zones in Figure 4c).Finally, the 3D edge enhancement method was used to enhance the edge detected in the variance attribute (Figure 4d).

Ant Tracking Result
In order to obtain good results, the ants should be restricted in a certain direction.Figure 5a shows the ant tracking attribute without using any azimuthal filter.Results show that most ant tracking attributes are continual along the W-E and N-S directions.However, in the study area, most faults are nearly W-E distributed, as shown in Figure 2. Hence, ant tracking was processed using two directional azimuthal filters, e.g., W-E and N-S directions with an angle of 60° as the

Ant Tracking Result
In order to obtain good results, the ants should be restricted in a certain direction.Figure 5a shows the ant tracking attribute without using any azimuthal filter.Results show that most ant tracking attributes are continual along the W-E and N-S directions.However, in the study area, most faults are nearly W-E distributed, as shown in Figure 2. Hence, ant tracking was processed using two directional azimuthal filters, e.g., W-E and N-S directions with an angle of 60 • as the tolerance.Figure 5b,c show the ant tracking attribute by using a W-E and N-S directional azimuthal filters, respectively.Results show that some W-E directional fractures were tracked very well by using the W-E directional azimuthal filter (Figure 5b); while the N-S directional fractures were tracked very well by using the N-S directional azimuthal filter (Figure 5c).These two ant tracking attributes were merged into one attribute (Figure 5d) which was used to extract fractures in Section 3.3.
Energies 2017, 10, 183 8 of 19 tolerance.Figure 5b,c show the ant tracking attribute by using a W-E and N-S directional azimuthal filters, respectively.Results show that some W-E directional fractures were tracked very well by using the W-E directional azimuthal filter (Figure 5b); while the N-S directional fractures were tracked very well by using the N-S directional azimuthal filter (Figure 5c).These two ant tracking attributes were merged into one attribute (Figure 5d) which was used to extract fractures in Section 3.3.5b) that did not disappear in Figure 5c; the blue ellipses indicate the detected N-S discontinuities (appear in Figure 5c) that were not disappeared in Figure 5b.

Extraction of the Large-Scale Fractures
The extraction result was strongly impacted by the parameter inputted in the extraction process.Through comparing the results of different parameters, we found that the extracted results were sensitive to deviation from a plane which was explained in Section 2.3.3.Figure 6 shows the extracted fault patches (coloured planes) by using different values of deviation from a plane parameter.The bold black lines represent the discontinuities detected by ant tracking.Comparing the four figures in Figure 6, the patches in the red ellipse in Figure 6b matches the ant tracking result, in which the value of deviation from a plane (DFAP) equals 7.There is another feature that the bigger value of DFAP we use, the more extended the patches will be (pointed with red arrows in Figure 6).As the value of DFAP is enlarged to 11, the patch even extends to another patch and  5b) that did not disappear in Figure 5c; the blue ellipses indicate the detected N-S discontinuities (appear in Figure 5c) that were not disappeared in Figure 5b.

Extraction of the Large-Scale Fractures
The extraction result was strongly impacted by the parameter inputted in the extraction process.Through comparing the results of different parameters, we found that the extracted results were sensitive to deviation from a plane which was explained in Section 2.3.3.Figure 6 shows the extracted fault patches (coloured planes) by using different values of deviation from a plane parameter.The bold Energies 2017, 10, 183 9 of 19 black lines represent the discontinuities detected by ant tracking.Comparing the four figures in Figure 6, the patches in the red ellipse in Figure 6b matches the ant tracking result, in which the value of deviation from a plane (DFAP) equals 7.There is another feature that the bigger value of DFAP we use, the more extended the patches will be (pointed with red arrows in Figure 6).As the value of DFAP is enlarged to 11, the patch even extends to another patch and merges with it (showed in a black circle in Figure 6d).The number of the extracted fractures decreases as the value of DFAP increases, which prove that some fractures are merged when we use the bigger value of DFAP.Hence, from the analysis of Figure 6, we selected the DFAP value equals 7 in this study.
Energies 2017, 10, 183 9 of 19 merges with it (showed in a black circle in Figure 6d).The number of the extracted fractures decreases as the value of DFAP increases, which prove that some fractures are merged when we use the bigger value of DFAP.Hence, from the analysis of Figure 6, we selected the DFAP value equals 7 in this study.Figure 7 shows the histogram of the length of the extracted fractures.There are a total of 1264 fractures.Results show that the fracture frequency increases with a decreasing fracture length; most fracture lengths range from 150 m to 1000 m and the average fracture length is about 511 m.The cumulative probability distribution of fracture lengths was fitted by using an exponential equation.Results show that there is a good fitting for fractures with length longer than 200 m (indicated by an arrow in Figure 7); the distribution of fractures which shorter than 200 m significantly deviated from the exponential distribution; we supposed that the deviation is because the fractures shorter than 200 m is beyond the identification capability of ant track using the seismic data.Hence, 200 m was used as a threshold for the division of small-scale fractures (SSFs) through assuming the large-scale fractures (LSFs).Figure 7 shows the histogram of the length of the extracted fractures.There are a total of 1264 fractures.Results show that the fracture frequency increases with a decreasing fracture length; most fracture lengths range from 150 m to 1000 m and the average fracture length is about 511 m.The cumulative probability distribution of fracture lengths was fitted by using an exponential equation.Results show that there is a good fitting for fractures with length longer than 200 m (indicated by an arrow in Figure 7); the distribution of fractures which shorter than 200 m significantly deviated from the exponential distribution; we supposed that the deviation is because the fractures shorter than 200 m is beyond the identification capability of ant track using the seismic data.Hence, 200 m was used as a threshold for the division of small-scale fractures (SSFs) through assuming the large-scale fractures (LSFs).
Figure 8 shows the dip angle and azimuthal distribution of the LSFs.Results show that the dip angle for most LSFs ranges from 0 • to 90 • with an average of 44.5 • .The dip azimuth can be grouped into the NW-SE direction.The distribution of the azimuth can be used as the parameter to generate the SSFs by assuming the LSFs and the SSFs follow the same characteristic of the azimuth distribution.
Figure 9 shows the histogram of the fracture dip angle which was measured from the 223 m core of seven wells.Results show that core-observed fractures mostly show a dip angle ranging from 40 • to 90 • .The average dip angle is 64 • .
Results show that there is a good fitting for fractures with length longer than 200 m (indicated by an arrow in Figure 7); the distribution of fractures which shorter than 200 m significantly deviated from the exponential distribution; we supposed that the deviation is because the fractures shorter than 200 m is beyond the identification capability of ant track using the seismic data.Hence, 200 m was used as a threshold for the division of small-scale fractures (SSFs) through assuming the large-scale fractures (LSFs).

Core Observation Result
Figure 10 shows the comparison of wireline logs, lithology, and fracture density at well A-67, whose location is shown in Figure 2. The fracture density was calculated by using the total area of the fractures on the core divided by the volume of the core (P32).In the column of "FracDen", the observed fracture density is presented using different colours.

Core Observation Result
Figure 10 shows the comparison of wireline logs, lithology, and fracture density at well A-67, whose location is shown in Figure 2. The fracture density was calculated by using the total area of the fractures on the core divided by the volume of the core (P32).In the column of "FracDen", the observed fracture density is presented using different colours.

Fracture Density Distribution
Table 4 lists the core observed fracture density and the corresponding values for AC, DEN, CN, GR, and ant attribute.The ANN model, which was explained in Section 2.3.4,was used to predict the fracture density at boreholes.Figure 11 shows the comparison of observed fracture density and the predicted fracture density.The correlation coefficient is about 0.60 which indicates a good match between them.

Fracture Density Distribution
Table 4 lists the core observed fracture density and the corresponding values for AC, DEN, CN, GR, and ant attribute.The ANN model, which was explained in Section 2.3.4,was used to predict the fracture density at boreholes.Figure 11 shows the comparison of observed fracture density and the predicted fracture density.The correlation coefficient is about 0.60 which indicates a good match between them.Figure 12 shows the distribution of AC, DEN, CN, GR, and ant tracking attribute, which were used to predict the fracture density distribution in 3D by using the built ANN model.Figure 13 shows the fracture density distribution.Results show that fracture density tends to be larger at the top of the erosional surface and near the faults.It is worth noting that faults are not necessarily flow channels or flow barriers because of their different sealing abilities; areas of more intense fractures are likely to be more conductive than areas of less intense fractures.Figure 12 shows the distribution of AC, DEN, CN, GR, and ant tracking attribute, which were used to predict the fracture density distribution in 3D by using the built ANN model.Figure 13 shows the fracture density distribution.Results show that fracture density tends to be larger at the top of the erosional surface and near the faults.It is worth noting that faults are not necessarily flow channels or flow barriers because of their different sealing abilities; areas of more intense fractures are likely to be more conductive than areas of less intense fractures.

DFN Modelling
Deterministic and stochastic methods were used in this study, respectively.The deterministic fracture model was generated by converting the extracted fracture patches from seismic data to the DFN model, which are LSFs.Table 5 lists the statistical parameters for the LSFs of two different fracture sets that are divided by azimuth.Results show that the dip angle, concentration, and anisotropy are similar for the two fracture sets.The stochastic fracture model was generated based on the spatially distributed fracture density, which was explained in Section 3.5, statistical data of fracture length, the maximum length of implicit fractures (50 m), and fracture orientation.The statistical parameters of fracture length were calculated by assuming the SSFs follows the same distribution of fracture length of LSFs as shown in Figure 7.The fractures were shaped in rectangles with the elongation ratio (the ratio of the horizontal length to its vertical length) of 2.9.
Figure 14a shows the 3D deterministic model of LSFs, in which the fractures are displayed using different colours for their length.For the SSF model, a small area was shown in Figure 14b because it takes a long time to generate and display the SSF model for the whole area.Figure 14c shows the fracture density distribution, which was used to generate the SSF model, for the same area as Figure 14b.
Malin et al. [60] reported that fracture-governed fluid flow in crustal rock is spatially correlated at all scale lengths.In this study, we did not assume any correlations between fractures because (1) we do not have enough information to study the correlation; and (2) the correlation studying is beyond the scope of this study.In the future, the spatially correlation of these large-and small-scale fractures will be studied.

DFN Modelling
Deterministic and stochastic methods were used in this study, respectively.The deterministic fracture model was generated by converting the extracted fracture patches from seismic data to the DFN model, which are LSFs.Table 5 lists the statistical parameters for the LSFs of two different fracture sets that are divided by azimuth.Results show that the dip angle, concentration, and anisotropy are similar for the two fracture sets.The stochastic fracture model was generated based on the spatially distributed fracture density, which was explained in Section 3.5, statistical data of fracture length, the maximum length of implicit fractures (50 m), and fracture orientation.The statistical parameters of fracture length were calculated by assuming the SSFs follows the same distribution of fracture length of LSFs as shown in Figure 7.The fractures were shaped in rectangles with the elongation ratio (the ratio of the horizontal length to its vertical length) of 2.9.
Figure 14a shows the 3D deterministic model of LSFs, in which the fractures are displayed using different colours for their length.For the SSF model, a small area was shown in Figure 14b because it takes a long time to generate and display the SSF model for the whole area.Figure 14c shows the fracture density distribution, which was used to generate the SSF model, for the same area as Figure 14b.
Malin et al. [60] reported that fracture-governed fluid flow in crustal rock is spatially correlated at all scale lengths.In this study, we did not assume any correlations between fractures because (1) we do not have enough information to study the correlation; and (2) the correlation studying is beyond the scope of this study.In the future, the spatially correlation of these large-and small-scale fractures will be studied.

Conclusions
An efficient workflow of the discrete fracture network modelling in a buried hill reservoir was presented in this study.The results show that: (1) The ant tracking attribute is feasible to detect the small faults and large-scale fractures in the buried hill reservoir.Before using the ant tracking algorithm, the seismic condition processes are necessary to enhance the quality of the seismic data; (2) The variance attribute is a powerful method to detect the edge of the amplitude volume, which is essential to the ant tracking process.Particularly, we tried an automatic workflow to extract the faults and fractures detected by the ants.The automatical workflow can save computing time, although we lost some accuracy since the shape and the position of the fractures do not always follow the discontinuity of the ant tracking attribute; (3) In total, 1264 subtle faults and LSFs were extracted in the reservoir.These fractures can be divided into two sets in terms of the azimuth distribution.Mostly, the strike of the fractures is in the W-E direction; (4) The fracture density is predicted through an integrated method by combining the core observation, well log interpretation, and artificial neural network algorithm.The spatial distribution of the fractures matched the faults.

Conclusions
An efficient workflow of the discrete fracture network modelling in a buried hill reservoir was presented in this study.The results show that: (1) The ant tracking attribute is feasible to detect the small faults and large-scale fractures in the buried hill reservoir.Before using the ant tracking algorithm, the seismic condition processes are necessary to enhance the quality of the seismic data; (2) The variance attribute is a powerful method to detect the edge of the amplitude volume, which is essential to the ant tracking process.Particularly, we tried an automatic workflow to extract the faults and fractures detected by the ants.The automatical workflow can save computing time, although we lost some accuracy since the shape and the position of the fractures do not always follow the discontinuity of the ant tracking attribute; (3) In total, 1264 subtle faults and LSFs were extracted in the reservoir.These fractures can be divided into two sets in terms of the azimuth distribution.Mostly, the strike of the fractures is in the W-E direction; (4) The fracture density is predicted through an integrated method by combining the core observation, well log interpretation, and artificial neural network algorithm.The spatial distribution of the fractures matched the faults.

Figure 1 .
Figure 1.The workflow for DFN modelling used in this study.Orange areas represent original data; light blue areas represent processed data or methods; the red area represents results.

Figure 1 .
Figure 1.The workflow for DFN modelling used in this study.Orange areas represent original data; light blue areas represent processed data or methods; the red area represents results.

Figure 2 .
Figure 2. Locations of the Jingbei Oilfield and wells (a), and stratigraphy (b), in the study area.

Figure 2 .
Figure 2. Locations of the Jingbei Oilfield and wells (a), and stratigraphy (b), in the study area.

Figure 3 .
Figure 3.A graphical illustration of the ANN model used in this study.

Figure 3 .
Figure 3.A graphical illustration of the ANN model used in this study.

Figure 4 .
Figure 4.The result of the seismic data condition.The original seismic amplitude volume (a); the amplitude volume after structure smoothing (b); the variances attribute volume calculated after structure smoothing (c); and the 3D edge enhancement result to enhance the edge of the variance volume (d).

Figure 4 .
Figure 4.The result of the seismic data condition.The original seismic amplitude volume (a); the amplitude volume after structure smoothing (b); the variances attribute volume calculated after structure smoothing (c); and the 3D edge enhancement result to enhance the edge of the variance volume (d).

Figure 5 .
Figure 5.A time slice showing the ant tracking attribute.The results of ant tracking without azimuthal filter (a); ant tracking using a N-S filter (b); ant tracking using a W-E filter (c); the ant tracking result was merged with the ant tracking results extracted using different azimuthal filters (d).The red ellipses indicate the detected W-E discontinuities (appear in Figure5b) that did not disappear in Figure5c; the blue ellipses indicate the detected N-S discontinuities (appear in Figure5c) that were not disappeared in Figure5b.

Figure 5 .
Figure 5.A time slice showing the ant tracking attribute.The results of ant tracking without azimuthal filter (a); ant tracking using a N-S filter (b); ant tracking using a W-E filter (c); the ant tracking result was merged with the ant tracking results extracted using different azimuthal filters (d).The red ellipses indicate the detected W-E discontinuities (appear in Figure5b) that did not disappear in Figure5c; the blue ellipses indicate the detected N-S discontinuities (appear in Figure5c) that were not disappeared in Figure5b.

Figure 6 .
Figure 6.The extraction result of the discontinuity from the ant tracking attribute.The value of deviation from plane = 5 (a); the value of deviation from plane = 7 (b); the value of deviation from plane = 9 (c); value of deviation from plane = 11 (d).

Figure 7 .
Figure 7. Distribution of the length of the extracted fractures.The distribution of fracture length corresponds the exponential law.The expression of the fitting curve is y = 155.22e−0.054x , where x indicates the fracture length and y indicates the cumulative probability density.

Figure 6 .
Figure 6.The extraction result of the discontinuity from the ant tracking attribute.The value of deviation from plane = 5 (a); the value of deviation from plane = 7 (b); the value of deviation from plane = 9 (c); value of deviation from plane = 11 (d).

Figure 7 .
Figure 7. Distribution of the length of the extracted fractures.The distribution of fracture length corresponds the exponential law.The expression of the fitting curve is y = 155.22e−0.054x , where x indicates the fracture length and y indicates the cumulative probability density.

Figure 7 . 19 Figure 8
Figure 7. Distribution of the length of the extracted fractures.The distribution of fracture length corresponds the exponential law.The expression of the fitting curve is y = 155.22e−0.054x , where x indicates the fracture length and y indicates the cumulative probability density.

Figure 8 .
Figure 8. Rose diagram showing the distribution of the azimuth of the extracted fractures from seismic activity.

Figure 9
Figure9shows the histogram of the fracture dip angle which was measured from the 223 m core of seven wells.Results show that core-observed fractures mostly show a dip angle ranging from 40° to 90°.The average dip angle is 64°.

Figure 9 .
Figure 9. Histogram of fracture dip angle from the 223 m core observation data from seven wells (a); a typical core photo showing fractures (b).

Figure 8 . 19 Figure 8
Figure 8. Rose diagram showing the distribution of the azimuth of the extracted fractures from seismic activity.

Figure 8 .
Figure 8. Rose diagram showing the distribution of the azimuth of the extracted fractures from seismic activity.

Figure 9
Figure9shows the histogram of the fracture dip angle which was measured from the 223 m core of seven wells.Results show that core-observed fractures mostly show a dip angle ranging from 40° to 90°.The average dip angle is 64°.

Figure 9 .
Figure 9. Histogram of fracture dip angle from the 223 m core observation data from seven wells (a); a typical core photo showing fractures (b).

Figure 9 .
Figure 9. Histogram of fracture dip angle from the 223 m core observation data from seven wells (a); a typical core photo showing fractures (b).

Figure 11 .
Figure 11.Observed fracture density against the predicted fracture density by the ANN model.

Figure 11 .
Figure 11.Observed fracture density against the predicted fracture density by the ANN model.

Figure 13 .
Figure 13.Vertical profile of fracture density.Black lines indicate the faults.The location section is shown in the small inset figure, which shows the fracture density distribution on the top grid in the model.

Figure 13 .
Figure 13.Vertical profile of fracture density.Black lines indicate the faults.The location section is shown in the small inset figure, which shows the fracture density distribution on the top grid in the model.

Figure 14 .
Figure 14.Discrete fracture model of large-scale fractures and small-scale fractures.The model of the LSFs (a); the red square indicates the area where SSFs were generated.The model of SSFs (b); the fracture density of the SSFs are in the same position as Figure 14b (c).

Figure 14 .
Figure 14.Discrete fracture model of large-scale fractures and small-scale fractures.The model of the LSFs (a); the red square indicates the area where SSFs were generated.The model of SSFs (b); the fracture density of the SSFs are in the same position as Figure 14b (c).

Table 1 .
List of parameters used in conditioning the seismic volume.

Table 2 .
The ant tracking parameters used in this study.

Table 3 .
Parameters used in the extraction of the fault patches in this study.

Table 3 .
Parameters used in the extraction of the fault patches in this study.

Table 4 .
Core observed fracture density and the corresponding wireline log values and ant attribute. AC,

Table 4 .
Core observed fracture density and the corresponding wireline log values and ant attribute.

Table 5 .
Statistical parameters of different fracture sets.

Table 5 .
Statistical parameters of different fracture sets.