Fluid flow is strongly affected by fractures in unconventional reservoirs. It is essential to deeply understand the flow characteristics with fractures for improving the production and efficiency of unconventional reservoir exploitation. The purpose of this work is to develop an accurate numerical model to evaluate the transient-pressure response for well intersecting fractures. The meshes generated from Fullbore Formation Micro-Imager (FMI) images ensure an efficient numerical description of the geometries for fractures and interlayers. The numerical simulation is implemented by an inhouse finite element method-based code and benchmarked with drill stem test (DST) data. The results show that three flow regimes appear in the reservoir with fractures within the test period: wellbore afterflow, pseudolinear flow, and radial flow. In contrast, only the wellbore afterflow and radial flow appear for the wells without fractures. The results also reveal that fractures dominate the flow near the wellbore. Verification and application of the model show the practicability of the integrated approach for investigating the transient-pressure behaviors in the unconventional reservoir.

Many fractures have formed in coal after the long-term geological process, which significantly influences coal seam gas (CSG) development [1, 2]. Multiple scale fractures have been investigated, ranging from nanoscale to field scale with hundreds of meters [3–6]. Microscale fracture analysis is mainly based on focused ion beam scanning electron microscopy (FIB-SEM) [7], scanning electron microscopy (SEM) [8], nuclear magnetic resonance (NMR) [9], X-ray micro-computed tomography (Micro-CT) [10]. Field-scale fractures can be detected by geophysical methods [11]. However, the characterization of meter-scale fractures in coal seams directly is still a challenge.

In petroleum engineering, the transient-pressure response is usually used to estimate the flow resistance, fault, and boundary effect of the formation, which can also be used to characterize fractures [12]. A large number of analytical models for wells with vertical fractures and horizontal fractures have been proposed [13–16]. However, the fracture pattern information is arbitrary in most cases [17]. In recent years, the analytical model of wells with inclined fractures has appeared gradually. Dinh and Tiab [18] developed a fully penetrating inclined fracture model based on uniform flux and infinite conductivity. Teng and Li [19] divided fracture into small panels, each used as a panel source. On this basis, they developed a semianalytical model of a finite-conductivity partially penetrating inclined fracture and distinguished the flow regimes. Generally speaking, the commonly used analytical models are overly simplified and idealized. All models are based on a single fracture or predefined fracture pattern or fracture network, which makes it impossible to reflect the complex reservoirs in reality [17]. Due to the natural limitations of analytical or semianalytical models, the numerical models, which can deal with boundary and fault more conveniently and simulate inclined well and horizontal well flexibly, are gradually developed [20]. In the past decades, a large number of numerical models have been established to study the transient pressure response and flow characteristics [21–25]. Along this line, Kamal et al. [26] proposed a systematic method named Numerical Well Testing (NWT) consisting of five steps: conventional analysis, locally refining grids, modifying numerical models, upscaling, and validating the models. Matthäi et al. [27] developed a fluid-flow model considering the geometry, thickness, and layered permeability of the normal faults for transient well testing. Agada et al. [28] established a high-resolution three-dimensional (3D) outcrop geological model, which was used for well test analysis and flow simulation after upscaling and fracture modeling. Yan et al. [29] presented an embedded discrete fracture model based on mimetic finite difference method, which can be calculated with permeability tensor and simulated complex fractured reservoir. Pan et al. [30] proposed a systematic workflow to integrate well testing data of a single well in a discrete fracture/matrix (DFM) model or a dual-porosity dual-permeability (DPDK) model for optimizing recovery. Pouladi et al. [31] developed a numerical model considering complex fault structures to investigate fluid flow characteristics. However, the above models do not fully view the actual and specific fracture information near the wellbore.

With the development of computer imaging technology, fractures and interlayers in reservoirs can be characterized accurately through the FMI log. In this work, we have developed a numerical model based on the FMI image of well Lake Goran 1 in Gunnedah basin to investigate the pressure dynamics of CSG reservoirs. We compare the simulation data with DST data to verify the accuracy of the model. Then, with the aid of the proposed numerical model, we analyze the influence of different fracture parameters on the flow regimes. Subsequently, we use the model to investigate the flow and pressure evolution of fractures and matrices. The novelty of this study is the fine characterization of fractures in geological model by FMI image and the organic unity of well logging and well testing.

The Gunnedah Basin, New South Wales, Australia, is a Cisuralian to Middle Triassic structural basin covering an area of over 115,000 km2 [32–34]. It contains up to 1200 meters of Permian and Triassic sediments that hosts tremendous thermal coal and CSG resources [35]. The basin has more than 500 Gt of coal and considerable potential reserves of coalbed methane [32]. The Gunnedah basin is unconformably located on the Permian or older basement and partially unconformably covered by sedimentary and volcanic rocks of the Jurassic and Cretaceous Surat Basin sequence [36]. The large-scale magmatic intrusion occurred in Late Triassic to Early Jurassic [37, 38]. Figure 1 shows the geological map of the basin and the position of the target well used for numerical simulation.

The study borehole, Lake Goran 1, is located in the southeast of the basin. The total depth of the well is 945.7 m. The well passes through Black Jack Group, Millie Group, and Bellata Group from top to bottom. The CSG targets we are concerned about are in the Hoskissons Coal in the Black Jack Group (405.4 m–414.1 m; Figure 2). The target horizon consists of three interlayers, eight open fractures, and two closed fractures [34]. Logging data, DST data, and core photos can be obtained from the DIGS website, a public online archive maintained by the Geological Survey of NSW (GSNSW).

This study aims to generate a formation model by using FMI image and then investigate the influence of fractures on well test curve and fluid flow characteristics. Although a considerable number of analytical and numerical models have been developed before, few models can make full use of the detailed and accurate information of interlaminar and cracks provided by FMI images. The unique feature of this model is that it can fully connect well logging information with well test and reservoir development. In addition, the model has good scalability. Images of different scales, from core scanning images in the laboratory to FMI images of wells, can be used to study fluid flow.

The steps of using FMI images to generate reservoirs are as follows: firstly, the Criminisi algorithm or deep learning algorithm is used to fill the blank strips of the image [39, 40]. If the image is complete, skip this step. Secondly, adjust the aspect ratio of the FMI image to the same scale as the mine. Thirdly, the inhouse Meshtool is used for generating meshes based on the image. Fourthly, keeping the vertical direction unchanged, the meshes are rolled into a wellbore shape along the radial direction. Fifthly, the duplicate nodes at the junction of the left and right meshes are deleted, and the remaining nodes are renumbered. Sixthly, the reservoir is formed by extending the meshes along the radial direction of the wellbore. It is worth noting that the meshes of the near-wellbore zone need to be highly dense, while the meshes far away from the wellbore can be sparse. This is because the pressure of the near-wellbore zone changes dramatically. Only high dense meshes can ensure the accuracy of numerical simulation. Seventhly, the elements of interlayers and fractures are determined by the information provided by FMI image. Eighthly, adjust the parameters of fracture, interlayers, and matrix until DST or well test data can be fitted. Lastly, the developed model is applied to numerical simulation.

Next, a fracture will be taken as an example to introduce how to determine the fracture elements in detail based on the FMI image. The inclination angle (⁠θ⁠) of a fracture can be defined using Eq. (1):
(1)θ=arctanh/d,
where h is the vertical distance between the highest point and the lowest point of the fracture in the FMI image (Figure 3(a)) and d is the wellbore diameter (Figure 3(b)).
The orientation of the lowest point of the fracture represents the azimuth (⁠β⁠) of the fracture (Figure 3(a)). The coordinates (⁠xf⁠, yf⁠, zf⁠) of any point on the fracture in a three-dimensional space can be obtained. Then, the normal unit vector of fracture surface can be expressed as (⁠sinβ∙sinθ,cosβ∙sinθ,cosθ⁠); the formula of the fracture plane can be defined by Eq. (2):
(2)a∙x+b∙y+c∙z−a∙xf+b∙yf+c∙zf=0,
where a=sinβ∙sinθ⁠, b=cosβ∙sinθ⁠, c=cosθ⁠.

The elements we built in the model are hexahedrons, so each element has eight nodes. The coordinates of 8 nodes are expressed as x1,y1,z1⁠, x2,y2,z2⁠, ⋯⁠, x8,y8,z8⁠. The coordinates of the eight nodes are substituted into the left part of Eq. (2). If the partial results are more significant than 0 and the other partial results are less than 0, the element is penetrated by fracture. If all the results are more significant than 0 or all are less than 0, it means that the element is not penetrated by fracture. Do the above calculation for each element in the model until all fracture elements are found. The schematic diagram of fracture elements identification is shown in Figure 4.

4.1. Mathematical Theories

The inhouse finite element method-based code PANDAS is an advanced multiphysical field coupling software and has been applied in various scenarios including the interacting fault system dynamics and geothermal and unconventional reservoir analysis [41–43]. The fluid module of PANDAS is extended and applied to investigate the transient-pressure behaviors and fluid flow characteristics in a CSG reservoir here [43].

PANDAS is used to solve partial differential equations, including continuity equations and motion equations. The continuity equation of fluid flow in porous media can be expressed as follows:
(3)ϕct∂P∂t+∇∙ρv=0,
where ϕ is porosity; ct is the total compressibility; P is the fluid pressure; ρ is the fluid density; and v is the fluid velocity. ct=ϕcl+1−ϕcr⁠, where cl and cr are the compressibility of fluid and rock, respectively. In this paper, bold symbols are used for vector, tensor, and matrix variables.
Darcy’s law can express the motion equation of fluid in porous media:
(4)v=−kμ∇P−ρg∇D,
where k is the permeability tensor of formation; μ is the dynamic viscosity; g is the gravity acceleration; and D is the depth.

4.2. Basic Assumptions

Some simplifications are made in this study at a reasonable computational cost. The model established in this paper is based on the following assumptions:

  • (1)

    Gravity effect is negligible because the thicknesses of the model are less than 10 m

  • (2)

    Adsorption and desorption characteristics are ignored in the CSG reservoir

  • (3)

    The permeability values of matrix and fracture in x-, y-, and z-directions are the same

  • (4)

    Permeability does not change with pressure

  • (5)

    The mineral does not react with the injected fluid

  • (6)

    The wellbore storage effect is considered in the model

4.3. Model Parameters

The meshes created using the method in Section 3 are shown in Figure 5, which was applied to study the flow characteristics in the CSG reservoir. The reservoir height was 8.7 m, and the radius of wellbore and reservoir were 0.048 m and 42.5 m, respectively. The model was discretized into 104,760 nodes and 98,600 hexahedral elements. The model had three units. The first unit was matrix. The second unit was three interlayers that were almost impermeable. The third unit was eight inclined fractures. The top-down indexes of the eight open fractures in Figure 2 were marked as 1,2,⋯,8⁠. Then, the fracture information is shown in Table 1.

The boundary conditions of the model are shown in Figure 6. Fluid was produced at a constant flux rate qw⁠. Due to the wellbore storage effect, the flux rate q at the bottom hole in the early stage is usually less than qw (Figure 6(a)). The reservoir is cylindrical with upper, lower, and outer impermeable boundaries (Figure 6(b)). The basic parameters used in numerical simulation are shown in Table 2.

4.4. Validation of the Numerical Model

The developed numerical model is to be verified by DST data of Lake Goran 1. The whole DST process lasted for 6 hours. The initial flow lasted 10 minutes to release mud pressure, and then the initial shut-in lasted 50 minutes. The second flow lasted 1 hour to remove well plugging, and the second shut-in lasted 4 hours for pressure build-up curve. Since the model we built did not need to eliminate the influence of drilling fluid, we only had a flow and shut-in cycle. The flow lasted 1 hour in the model; then, the pressure build-up curve was recorded for 4 hours. The comparison between DST data and simulation data is shown in Figure 7. Obviously, the DST data and numerical simulation results agree well with each other.

5.1. Transient-Pressure Behavior

The purpose of our model is closely combined with the field well test. The numerical well test is set to last 300 hours, which is much longer than the current well test time of CSG reservoir and saves computing resources. In this section, we keep the geometry of fracture network and basic parameters (see Table 2) unchanged and only change fracture permeability and porosity to investigate the transient-pressure response. The flux rate equals 0.05 m3/day, and the specific experiment parameters are shown in Table 3.

Figure 8 plots the curves of pressure responses for different fracture parameters with the same matrix and interlayers. As one can see from the figure, fracture parameters have a significant influence on the transient-pressure as well as derivatives during the whole period. In experiment 1, the fracture parameters are utterly consistent with that of the matrix, which means that no efficient conductive channel exists in the formation. We can only distinguish two flow regimes from Figure 8(a) clearly: wellbore afterflow and radial flow. During the flow period of up to 300 hours, no other flow regimes such as linear flow appeared in the formation. Correspondingly, when efficient conductive channels exist in the formation, three flow regimes can be distinguished from Figure 8(b), Figure 8(c), and Figure 8(d): (1) wellbore afterflow, (2) pseudolinear flow, and (3) radial flow. Figure 9 demonstrates the comprehensive comparison of pressures and derivatives for the above four cases. As shown in Figure 9, greater fracture permeability leads to smaller pressure consumed due to the minor flow resistance in the formation. With the increase of fracture conductivity, the influence of pseudolinear flow is also gradually increasing. Figure 10 plots the four instantaneous pressure distributions marked in Figure 8 near the well, further showing three flow states in the formation.

In experiment 3, the permeability of fractures 1 to 7 is 3.00E−16 m2⁠, while the permeability of fracture 8 is 3.00E−15 m2⁠. Other parameters not shown in this table are identical in the four numerical experiments.

5.2. Flow Characteristics in Reservoir

We use the numerical simulation results of experiment 2 to discuss the fluid flow characteristics in the reservoir. Figure 11 shows the fluid flow velocity distribution near the well with the flow time of 2.5 hours, 11.7 hours, 29.0 hours, and 141.5 hours. In the model we developed, fluids can flow into the well from fractures and matrices. But due to the smaller flow resistance of fractures, the velocities of fluids in fractures are much higher than those in the matrices. Figure 11 shows the fluid flow evolution in the formation. A high-velocity zone is gradually formed around the fractures compared with the matrix farther away from the well. As the flow continues, the influence of fractures gradually weakens, and radial flow slowly appears in the formation.

We choose four points to investigate the evolution of ΔP and velocity in the flow process. We set the indexes of the five points as P1, P2, P3, P4, and P5 at the section Z=8.3 m with the distance from the wellbore of 0.43 m, 2.05 m, 3.00 m, 6.41 m, and 13.69 m. P1 is located on the fracture, P2 is on the edge of the fracture, and P3 to P5 are in the matrices and gradually away from the wellbore. The results are shown in Figure 12. During the flow period, the ΔP of the fracture increases rapidly, followed by the edge of the fracture, and then the ΔP of the points in the matrices increases in turn, which is inversely proportional to the distance from the point to the wellbore. Accordingly, the flow velocity of fracture is much higher than that of matrices. The above analysis reveals that fractures dominate the flow near the wellbore.

Subsequently, we investigate the streamlines during the flow period. The streamlines near the well at different sections Z=1.3 m⁠, Z=2.9 m⁠, Z=6.6 m⁠, and Z=8.3 m with the flow time of 141.5 hours are plotted in Figure 13. As we can see from the figure, the streamline paradigms at each section are widely different caused by the fracture combination indicated by Figure 5 and Table 1, which means that the flow field is dominated by the number, azimuth, and location of fractures. Moreover, the dense streamlines and the high velocity at the fractures also illustrate that the flow field is deeply affected by fractures.

In this work, a numerical well testing model for wells intersecting fractures has been proposed using the finite element method. The meshes of 104,760 nodes and 98,600 hexahedral elements generated from FMI images are successfully applied to describe the highly complex CSG reservoir with eight open inclined fractures and three intercalations. Based on fine network structure, the finite element simulation is benchmarked with DST data and used to investigate the transient-pressure response and flow characteristics in the CSG reservoir and further evaluate the influence of fractures on the evolution of pressure and velocity in the reservoir. Numerical simulation results show only wellbore afterflow and radial flow appear in 300 hours without fractures. The complex network of inclined fractures in the CSG reservoir exhibits three flow regimes: wellbore afterflow, pseudolinear flow, and radial flow. Furthermore, fractures dominate the flow near the wellbore with much higher flow velocity than that of the matrices. The proposed model will be extended with non-Darcy effect in the future to investigate the flow characteristics in the CSG reservoir more realistically.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors declare that they have no conflicts of interest.

This research work is funded by the National Natural Science Foundation of China (No. 52074251 and No. 92058211), the Fundamental Research Funds for the Central Universities (No. 202012003), the open fund (Grant Number SKLGP2017Z001) from State Key Laboratory of Geohazard Prevention and Geoenvironment Protection, and Ministry of Education of China (No. B20048).