Rapid productivity prediction method for frac hits affected wells based on gas reservoir numerical simulation and probability method

2.1 Embedded discrete fracture modeling

Fractal fracture network model based on fractal theory, extended finite element model based on finite element method, displacement discontinuous numerical model based on boundary element method, and finite element model based on finite element method are the main numerical simulation methods of fracture network in shale gas reservoirs. Additionally, there is the unconventional fracture network model built on the foundation of the discrete fracture network model, which is based on the dual medium model and the continuity principle. This led to the development of the EDFM, which combines the benefits of dual media and discrete fracture models [31,32,33,34,35,36,37,38,39,40,41].

This approach establishes a shale gas reservoir numerical model based on the gas reservoir numerical model in order to characterize the adsorbed gas and free gas in the shale gas reservoir. The liquid can flow between the matrix and matrix, and also can flow between the matrix and the fractures. The model considers Knudsen diffusion, adsorption and desorption, stress sensitivity, and other mechanisms, and uses the embedded discrete fracture treatment method to characterize the fracture network and natural fractures into the gas reservoir numerical model to generate multi-scale discrete fracture values for shale gas reservoirs simulation model. The reservoir is discretized using Cartesian grids, extra grids are introduced at the fractures to characterize mass transfer across the media, and the conductivity between the grids is used to quantify the effect of fractures on fluid flow.

Conductivity adopts non-adjacent connection method:

where q is the volume flow of phase l between non-adjacent connected pairs; λ _l is the relative mobility of phase l; T _NNC is the conductivity coefficient between non-adjacent connected pairs; ΔP is the pressure difference between grids; k _NNC is NNC penetration rate; A _NNC is the area of NNC; and d _NNC is the NNC distance.

1. NNC type I: It is the connection between the fracture segment and the matrix grid through which it passes.

   (3) T f-m = 2 A f ( K → → ⋅ n → ) ⋅ n → d f-m ,

   (4) d f-m = ∫ V x n d V V ,

   where A _f is the area on one side of the fracture section; K is the permeability tensor; n is the normal vector of the fracture plane; d _f−m is the average distance from the matrix to the fracture; V is the volume of the fracture grid; and x _n is the distance between the matrix unit and the fracture unit.

2. NNC type Ⅱ: It is the connection between different fracture segments of a single fracture.

   (5) T seg = T 1 T 2 T 1 + T 2 ,

   (6) T 1 = k f A c d seg 1 , T 2 = k f A c d seg 2 ,

   where k _f is the fracture permeability; A _c is the area of the common surface of the two fracture segments; d _seg1 is the distance from the center of the fracture segment 1 to the common surface; d _seg2 is the distance from the center of the fracture segment 2 to the common surface.

3. NNC type Ⅲ: It is the connection between intersecting fractured sections.

where T _int is the conductivity between the intersecting fractures; L _int is the length of the intersection line; k _f1 and k _f2 are the permeability of fractures 1 and 2, respectively; w _f1 and w _f2 are the openings of fractures 1 and 2, respectively; d _f1 and d _f2 are, respectively, the weighted average distance from the sub-segment of fractures 1 and 2 to the intersection line of the fracture.

2.2 Automatic history matching of gas reservoir numerical models

The modeling process involves numerous, intricate, and unknown aspects, which cause variations between the gas reservoir model and the actual gas reservoir. The numerical model of the gas reservoir must be corrected and optimized based on the production data history matching, and the reservoir parameters are then reversed and corrected in accordance with the observed actual gas reservoir (well) performance. The automatic history matching technology of the gas reservoir numerical model overcomes the shortcomings of manual trial calculations, and uses the optimization method to automatically correct the model parameters and structure, reducing the fitting time and achieving higher accuracy.

Aiming at the established numerical model of shale gas reservoir, this method adopts the ES-MDA automatic history matching technology to carry out the history matching process of the model [42]. The algorithm updates the model parameters based on the difference between the simulated data and the historical observation data. The multi-iteration set smoothing algorithm is iteratively updated and solved according to the correlation between the simulated data. It is suitable for high-dimensional situations. With a small number of iterations, the multi-iteration set smoothing algorithm can obtain a better solution.

After completion of the production data history matching gas reservoir numerical model, it is possible to invert fracture conductivity, fracture half-length, reservoir pressure, reservoir permeability, and gas saturation in the artificial fracture network, based on production performance data and other important geological and engineering factors.

2.2.2 Autoencoder

The ES-MDA algorithm can iteratively update the solution for high-dimensional variables, but it involves the inverse process of model parameters and production data, and the high-dimensional variables will increase the computational consumption, so the model parameters are dimensionally reduced using autoencoder neural network.

Autoencoder neural network is an unsupervised machine learning method that enables data reconstruction or dimensionality reduction through encoding and decoding operations. Its neural network structure is shown in Figure 1. Among them, x represents the input sample data, h is the hidden layer node, and x' is the output data. The loss function of the autoencoder neural network can be expressed as the square of the norm of the difference between the input and output.

(13) ℒ ( x , x ' ) = ∥ x − x ' ∥ 2 .

Figure 1

Structure of the autoencoder.

Through the training of the neural network, the best effect is achieved in the test samples, which can realize the encoding process from the input to the hidden layer, i.e., data dimensionality reduction, and the decoding process from the hidden layer to the output, i.e., data reconstruction. The samples are randomly generated according to the range of shale gas reservoir model parameters, the autoencoder neural network is trained, the samples are input to the self-encoder to get the hidden layer data, and the hidden layer data are used as parameter variables to achieve the purpose of model parameter dimensionality reduction.

The overall history matching flowchart is shown in Figure 2. First, the shale gas reservoir parameters are characterized, including fractured fractures, natural fractures, relative permeability, stress sensitivity, matrix, etc. Sample data are randomly generated according to the parameter range, and parameter downscaling is performed using a autoencoder to obtain the hidden layer data, which is used as the solution variable. Model parameters are obtained by reconstructing the hidden layer data, and reservoir numerical simulation is performed to obtain the model simulated production data. Termination conditions are judged, and if the termination conditions are met, the history fitting ends; if the termination conditions are not met, the model parameters are updated using an ensemble smoothing algorithm until termination.

Figure 2

Workflow of history matching.

3.1 Probabilistic characterization of frac hits affected wells

Frac hits mainly occurs in the alternate layout of new and old wells. When a new well is being fractured, some of the fracturing fluid will affect nearby production wells, causing their output to drop sharply or cease entirely. The main factor contributing to pressure channeling is that the well that needs to be reformed in this area has been producing next to the old well for a while, which lowers the regional reservoir pressure and alters the reservoir’s original in situ stress field. When new wells undergo fracturing, due to changes in the regional reservoir pressure and in situ stress fields, artificial fracturing fractures are prone to propagate to the pressure-reducing area and cause fracturing fluid to interfere with the neighboring wells [43]. At the same time, part of the new and old wells’ natural fractures are developed, and when the new wells are reformed, they act as channels with high conductivity, which can cause fracturing fluids to flow to adjacent producing wells [44].

Due to uncertainty and the lack of data, the probability approach, also known as the uncertainty method, is frequently used to determine the probability distribution of oil and gas geological reserves. The probability distribution of recoverable reserves can also be computed using this technique in conjunction with the production decline approach and other dynamic methods. The probability technique uses Monte Carlo simulation to complete the probability distribution of the value of the objective function, considering the value of the probability distribution of the uncertainty parameter in the problem to be addressed.

After the occurrence of frac hits, the production of a single well is affected by fracturing fluid interference, and the production cannot be restored within a certain period of time, and its final recoverable reserves will also be different from before. During the period when the productivity of a single well is not restored, the prediction of the recoverable reserves of a single well becomes a probability distribution problem due to the influence of pressure channeling. During the process of frac hits between wells, the fracturing fluid flow to adjacent production wells through artificial fracturing fractures and natural fractures is disturbed. The water content in the wellbore of the disrupted production well and the nearby formation fractures increases significantly, resulting in water lock in the reservoir near the production well and fluid accumulation in the wellbore, which greatly affects the productivity of production wells. Fluid accumulation in the wellbore will cause changes in bottom hole flowing pressure, and reservoir damage will cause changes in the conductivity of the fracture network, regional permeability, and the number of fracture openings. The influence of pressure channeling on the uncertainty of the bottom hole flow pressure, fracture network conductivity, regional permeability, and the number of fracture openings turns the final recoverable reserves of a single well into a probability prediction problem.

Figure 3

Probability-based productivity prediction process for frac hits affected wells.

3.2 Method and process

1. Because of the uncertainty of the fracture network conductivity, the permeability of the reformed area, bottom hole flowing pressure, and the number of fracture openings in the process of pressure channeling first establish artificial fractures and natural fractures based on EDFM to characterize single wells. For matrix and micro-fractures in shale gas reservoirs, Multiple interacting continua model is used to characterize the unsteady crossflow between them. For hydraulic fractures formed by artificial fracturing, EDFM is used to explicitly characterize them. The explicit large fracture is directly embedded into the equivalent micro-fracture grid, and the large fracture elements are divided by the intersection point between the large fracture and the grid boundary. By constructing the quasi-steady crossflow between each large fracture element and the corresponding micro-fracture grid, the flow field coupling is realized.

2. The gas reservoir numerical model is then used to perform automatic history matching on the historical production data before frac hits based on the dimensionality reduction and ensemble smoother algorithm which is introduced in Section 2.2. Then, we can obtain the conductivity of the fracture network before frac hits, the permeability of the reformed area, the bottom hole flowing pressure, and the number of fracture openings.

3. The probability distribution function is used to quantify the uncertainty of these random variables with a certain range of values. The theoretical distribution model has specific functional expressions, including uniform distribution, triangular distribution, normal distribution, lognormal distribution, and so on. When the parameters can only obtain a range of values and are uniformly distributed in this range, the uniform distribution model is usually satisfied. When parameters have not only range values, but also most probable values, they can be described by triangular distribution. Normal distribution and lognormal distribution are common in the field of oil and gas reserves estimation, which are usually used to describe reservoir area, thickness, porosity, saturation, and other parameters. In this work, the conductivity of the fracture network, the permeability of the reformed area, the bottom hole flow pressure, and the number of fractures opened are quantified using a triangular distribution function, and the permeability of the reservoir outside the reformed area is quantified using a normal distribution function.

4. In the end, based on the quantization function of the uncertainty parameter after frac hits, Monte Carlo simulation is used to extract the probability distribution of the uncertainty parameter in the single well production decline model, and bring it into the model for calculation. After completing multiple iterations, the predicted probability distribution of the final recoverable reserves is obtained, the probability values of P10, P50, and P90 are read from it, and the result of P50 is taken as the final recoverable reserves of the current single well (Figure 3).

4.1 Establishment and matching of gas reservoir numerical model

A regional gas reservoir numerical simulation model of Well W5 was created using the embedded discrete model modeling technique, and a geological model was developed utilizing geologically related knowledge results (Figure 5). The model covers the area where Well W5 is located in the lateral direction of 3,000 m × 800 m, the length of the coarse grid is 50 m × 50 m in the X and Y directions, and the average thickness in the longitudinal direction is 45 m. It is divided into seven layers according to the thickness of small layers. Under the given production conditions of gas production, the changes in geological engineering parameters of 215 days of well-opening production are simulated. Through the multiphase flow theory in the wellbore, after converting the historical production data of casing pressure into bottom hole pressure, the automatic history matching technology of gas numerical reservoir simulation is used to fit the production data and simulated data of bottom hole pressure (Figure 6). The parameters set during the matching process include fracturing fracture conductivity, matrix permeability, phase permeability curve, stress sensitivity curve, and other parameters, a total of 368 parameters. Set the number of iterations to 2, the set size to 100, and the total number of simulations is 200. After automatic history fitting, the error between the numerical model bottom hole pressure simulation data and the historical data gradually decreases, and finally reaches a stable state (Figure 7).

Figure 5

Pressure distribution before frac hits.

Figure 6

Matching curve of bottom hole pressure.

Figure 7

The error distribution of the objective function of automatic history matching.

After the gas reservoir numerical model is matched with the production data history, the relevant parameters are processed, and the fracture conductivity coefficient of the region near Well W5 before the occurrence of fracturing is 6.25 D*cm, the average equivalent permeability of the reformed zone is 0.0058 mD and the regional. The average permeability of the reservoir is 0.00023 mD. The bottom hole flowing pressure, fracture network conductivity coefficient, the permeability of the reformed area, and the number of fracture openings are brought into the triangular distribution function, and the permeability of the peripheral reservoir of the reformed area is brought into the normal distribution function. Based on the production data before frac hits, using the Monte Carlo simulation random sampling method and substituting the production decline model for 5,000 times, the probability distribution of the final recoverable reserves of Well W5 is obtained, and the final recoverable reserves of P50 are predicted to be 134 million cubic meters (Figure 8). Well W5 resumed production after 71 days of frac hit, and part of the production was restored. As a verification method, after resuming production for a period of time, using the Duong method (a = 6.2649, m = 1.4659), YM-SPED method (n = 0.2634, T = 1.4221), and WK method (λ = 0.0872) of the empirical production decline method, the single well EUR is calculated to be 141 million cubic meters, 133 million cubic meters, and 139.5 million cubic meters, respectively (Figure 9). The results are consistent with this article. The error of the prediction results of the method is 4.96, 0.75, and 3.94%, respectively. The accuracy of the prediction of this method has been well verified.

Figure 8

Probability method capacity prediction results.

Figure 9

Calculation results of three empirical production decline models. (a) Duong method predicted results. (b) YM-SPED method predicted results. (c) WK method predicted results.