This paper presents reservoir simulation modeling of a hydrocarbon accumulation with a fractured porous reservoir, incorporating the geomechanical effects of fracture closure during variations in formation pressure. The fracture permeability parameter is derived from the impact of stress on fracture walls. The fracturing parameter is determined based on 3D seismic data analysis. A permeability reduction model is implemented in the tNavigator reservoir simulation platform. The proposed approach improves the convergence of formation pressure dynamics in well data while maintaining flow rate and water cut adaptation accuracy. This results in enhanced formation pressure prediction and optimization of the pressure maintenance system.
During the development of a fractured reservoir, significant decreases in formation pressure trigger processes that are difficult to capture using single-medium models. Widely discussed in the literature, these processes include fracture closure and variations in fracture permeability, both of which affect the reservoir’s geomechanical state [1-3].
Currently, reservoir simulation models of oil and gas accumulations, even those containing fractured porous reservoirs, are often built using single-porosity models. These models are well-established and relatively straightforward to adapt to initial development conditions. However, well flow rates and bottomhole pressures depend heavily on the permeability and porosity of the near-wellbore formation. High water production is a function of phase permeabilities and aquifer activity. One of the most complex and crucial aspects of model calibration is the accurate simulation of formation pressure dynamics, as it impacts well performance, equipment selection, and stimulation planning. The well operating mode, the choice of downhole equipment, and the planning of inflow stimulation measures depend on formation pressure. From the point of view of reservoir physics, the content of free gas in the pores and, as a result, phase permeabilities, as well as fractures opening in the reservoir and absolute permeability depend on formation pressure.
The most accurate approach for modeling fractured reservoirs involves developing a discrete fracture network (DFN) model coupled with a dual-porosity/dual-permeability reservoir simulation. However, such models require extensive parametric input and are challenging to implement [4].
The article proposes an approach to geological reservoir simulation modelling of hydrocarbon accumulations with a fractured porous reservoir based on the geomechanical effects of the fracture space closure during changes in formation pressure. The permeability distribution along the section is dynamic. In the initial formation conditions, the reservoir is represented by both a pore medium and open fracture systems. With a decrease in formation pressure, an intensive decrease in fracture permeability occurs, as a result of which entire areas of the formation, previously considered a reservoir, cease to filter hydrocarbons. The initial value of permeability, as well as its change, are a function of stresses.
The study was conducted on the Tournaisian–Famennian formation (TFm) in one of the fields in the Perm Region. This formation consists of a massive carbonate bed with up to 19 layers and a total thickness of 50-75 m. The wells initially exhibit high flow rates (up to 100 m³/day or more), but most experience a several-fold decline in flow rates within the first 5-10 years.
To consider the dynamic behaviour of the permeability and porosity properties of a fractured reservoir during reservoir simulation modelling, modelling can be divided into several stages in accordance with the new approach:
To identify the effect of changing the filtration properties of the reservoir, the results of well flow studies as well as industrial data on well operation were analyzed at the first stage. The fact of a change in filtration properties can be preliminarily observed by a decrease in well permeability, determined by pressure recovery curves, by well indicator diagrams concave to the flow rate axis, as well as by the negative dynamics of the productivity coefficient [5, 6]. The first method of identifying this effect is the most reliable, since the permeability dynamics are analyzed on the well contour, where the colmation is absent. According to widespread practice, the dependence of permeability on pressure is approximated by an exponential function
where k0 and k – the initial and current permeability; p0 and p – the initial and current formation pressures; β – the permeability reduction degree.
The unknown values are the initial permeability under formation conditions k0 and the permeability reduction parameter β. Such an analysis is possible in the presence of a sufficiently wide range of formation pressures of well studies in non-stationary modes, at low water cut values and formation pressures exceeding the saturation pressure. At the pay under consideration, 10 wells correspond to these conditions, for which β was subsequently determined.
When selecting geological and geophysical criteria for the 3D distribution of the permeability reduction parameter, several approaches can be used. For example, in [7], ranking was performed by the P-wave travel time. The compressibility of the fracture walls was determined by the strength properties of the rock, which, in turn, depended on the P-wave velocity. The advantage of this method is the prevalence of acoustic logging materials for most pays. Its main disadvantage is the lack of reference to fracturing. In the case of a porous fractured reservoir, where fracturing is not widespread, based on the interval travel time of the P-wave, high rock compressibility can be falsely identified even in areas with a porous reservoir.
In this regard, the present work investigated the formation fracturing distribution based on well data and core studies. Since the number of wells covered by the fracturing studies is insufficient for simple interpolation, a search for a relationship between fracturing and 3D seismic materials was performed. A set of research data related to the identification of fractures was collected and analyzed [8]. Fracturing of eight wells in the TFm pay of the analyzed field was studied. The studies included a geophysical investigation of wells using the EMS by KarSar MS-110 and MicroScope devices, as well as fracturing determination in an oriented core.
The distribution of near-vertical fracturing with predominant strike azimuths of 90-130° and dip angles of 70-85° stands out from the entire data set. The results did not depend on the device used. The obtained distributions were further used to build a volumetric probabilistic model of fracturing in the studied pay. Classification was performed by machine learning methods in Python. The angle stack cubes of 3D seismic survey amplitudes were used directly (without conversion to attributes) as the initial parameters for training the model. The idea was to take a segment of the seismic trace passing through the studied point in the interval of ±50 m relative to this point using compression methods (principal vector methods, singular value decomposition, etc.) as input parameters of machine learning methods. Various boostings were used as machine learning models.
The classification quality was assessed by the cross-validation method using the ROC AUC metric and reached values of 0.88-0.91 for different models. Two classes (types) of void space were identified based on actual data: without fractures (0) and with fractures (1). Figure 1 shows the probability density distributions of the occurrence of these classes depending on the probability of predictions of machine learning models. The classes are well separated in such a way that with an increase in the predicted probability, the ratio of class 1 to class 0 increases.
Fig.1. Distribution of probability densities of fracturing for the selected classes (a) and the normalized difference in probability densities (b)
Fig.2. Probability of fracturing distribution along the roof of the TFm formation (a) and along section 1-1 (b)
Based on the obtained cube of a 3D fracturing probability distribution (Fig.2), the average value of this parameter in the perforation interval for all wells was determined. Further analysis of the parameter for a group of wells with a sufficient number of pressure build-up curves demonstrated its linear relationship with a decrease in permeability:
where CD – the probability parameter of the presence of fracturing (variation range from 0 to 1).
This dependence was subsequently used to identify regions with decreased permeability in the reservoir simulation modeling, as well as to compile tables of the dependence of permeability and porosity properties on formation pressure in each identified region.
At present, to obtain a 3D distribution of permeability in a reservoir simulation model of fractured and fractured porous pays, in most cases, correlations on porosity are still used with subsequent calibration to the results of well flow testing (WFT). A geomechanical approach was used to determine the permeability distribution.
There is a lot of information in scientific and technical literature about the influence of the stress state of rocks on the permeability of fractured porous and fractured reservoirs [9-11]. However, the permeability of porous reservoirs does not depend on the stress state as much as the permeability of fractured rock [12-14]. With this in mind, a search was performed for the relationship between the stress state of the massif [15, 16] and fracture permeability determined based on WFT.
It is clear that the permeability of fracture and fault systems is highly dependent on their orientation relative to the principal normal stresses. The critically stressed fault hypothesis [16, 17] states that “mechanically active faults are also hydraulically active, whereas mechanically inactive faults are not”, a view supported by numerous studies [16, 18, 19]. Mechanically active faults are under stresses that promote the opening of fractures, facilitating fluid flow. Thus, the permeability of a fracture system is influenced both by the normal stresses that compress it and by the tangential stresses that induce sliding [20]. Equations for estimating tangential and normal stresses acting along a fracture are widely known [1, 15].
The tendency for slip depends not only on these stresses but also on the rock’s strength. In general, the stress state criterion for a weakening surface (fracture or fault) is given by:
where μ – the friction coefficient; τ, σn – the tangential and normal stresses on the fracture surface, respectively.
Thus, having estimated the normal and shear stresses acting in the analysed fracture system and having compared them with the pay permeability data recorded in a certain interval, it is possible to judge the degree of the fracture system activity.
In the simplest model – assuming an impermeable matrix and a system of plane-parallel fractures – the fracture permeability is given by:
where e – the fracture aperture; d – the spacing between fractures.
As a parameter characterizing the fracture density (the reciprocal of the distance between fractures), the third degree of the fracturing probability can be used.
If we take the fracturing probability as a measure of fracture density (p = 1/d) and assume that the fracture aperture is proportional to the stress ratio (e = τ/σn), then an expression for the fractured reservoir permeability can be formulated as:
where Prob – the fracturing probability value; τ/σn – the ratio of normal (σn) and tangential (τ) stresses along the fracture system; a – an empirical coefficient [20].
In simulating fracture permeability, two types of fracturing were considered: ordered, corresponding to the maximum fracture concentration (azimuth 115°, dip angle 80°), and conditionally ordered, associated with the stratigraphic layer curvature, the strike azimuth of which is normal to the direction of maximum curvature. The maximum curvature azimuth was determined as a physical direction (rather than by comparing vectors along mutually perpendicular axes, as many geological modeling packages do). The dip angle for the conditionally ordered system was obtained by solving an optimization problem. Finally, the pore permeability from the geological model was added to the fracture permeability to yield the total reservoir permeability:
where Kpor – the pore component of permeability.
For subsequent flow simulations, the fracture component of permeability obtained in the local fracture coordinate system must be transformed into the global coordinate system of the model (2). For each fracture system, the permeability tensor is estimated in the global coordinate system of the reservoir simulation model [15]:
where α – the dip angle of a given fracture system; θ – the strike azimuth; m is the number of fracture systems.
The tangential and normal stresses and the directions in which they act in each cell for permeability simulation were derived from our geomechanical model. Developing a geological-geomechanical model [1, 18, 19] for a pay or an entire field involves assessing the initial stress state by determining parameters such as: vertical stress; minimum and maximum horizontal stresses; azimuth of horizontal stress; pore pressure; spatial distribution of the mechanical properties of both the pay and the surrounding rocks [21-23].
A 3D geomechanical finite element mesh was constructed by exporting the geometry (as eight-node element coordinates) from geological modeling software and converting it into a format compatible with the geomechanical simulator.
The grid cells had lateral dimensions of 100 × 100 m. In the region of the studied pay, the mesh was vertically refined to a cell thickness of 0.4 m. Instead of modeling the overlying strata in full, a leveling layer was created and loaded with a weight corresponding to the overburden during geomechanical computations.
Since the vertical stress is determined by the weight of the overlying rocks, the density distribution was simulated by interpolating density logging data and integrating over cell thicknesses to compute vertical stress. The mechanical properties of the pay were assigned based on relationships between static and dynamic characteristics obtained from core tests [24-27].
Lateral stress components were simulated using an isotropic elastic porous medium model [16, 28, 29]:
where Sh, SH – the minimum and maximum horizontal stresses, respectively; p – the formation (pore) pressure; σv – the vertical stress, Biot – the Biot coefficient; v – the Poisson ratio, E – the elastic modulus, ɛh, ɛH – the components of additional tectonic deformations.
The formation pressure distribution was adopted from the current reservoir simulation model and later adjusted using up-to-date pressure data (from well flow tests, hydraulic fracturing, etc.) [30].
Calibration of the simulated minimum horizontal stress profile was performed by comparing the estimated stresses with fracture closure pressures measured during hydraulic fracturing (Fig.3).
Subsequently, the finite element method (implemented in Ansys) was used to evaluate the stress-strain state of the rock mass. The simulation employed a coupled formulation that accounted for fluid flow in a deformable porous medium, with kinematic boundary conditions on the lower and lateral faces to incorporate tectonic deformations.
After completing the geomechanical computations (using equations (2)-(4), a permeability cube was generated and used in the reservoir simulation. The permeability values obtained by this method showed good agreement with WFT results (Fig.4).
Fig.3. Comparison of estimated stresses with measured values of hydraulic fracture closure pressures
1 – values in wells; 2 – line of equal values; 3 – regressive straight line
Fig.4. Crossplot of the cube root of the estimated and determined from WFT permeabilities
See Fig.3 for legend
In the final phase, reservoir simulation was conducted using the new permeability cubes and incorporating fracture closure concepts.
In the initial simulation, the reservoir was subdivided into many layers with limited hydrodynamic connectivity between interlayers. The fracture system-evident from core studies and microimager logs-was not explicitly modeled. Although some sections might formally be classified as nonreservoir (owing to negligible total porosity), fluid flow can occur via fractures.
To account for fracturing, we propose an intermediate model between the original single-medium model and a dual-medium (dual-porosity/dual-permeability) model with the following features:
Permeability reduction was implemented using 10 rocktab tables describing permeability change [31]. The regions of applying these tables were determined by the fracturing probability value. The first region corresponds to the accumulation areas with a minimum number of fractures, the tenth region – to the areas with the highest fracture density and maximum permeability reduction. The parameter β is determined by the linear dependence (1). The method for considering the effect of permeability reduction in flow simulations is implemented in the tNavigator software package.
Flow simulations showed that when the formation pressure decreases, there are areas where permeability drops especially rapidly. In the case of low initial permeability and a rapid decrease in permeability, zones may arise where fluid filtration ceases. This indicates an inconstant (dynamic) character of the lithological section over time (Fig.5).
The new approach improved the match between simulated and actual formation pressures in the reservoir model (Fig.6). For example, in well 3 the decline in formation pressure was accompanied by a drop in liquid flow rate from 80 to 20 m³/day. In the original model, a reduction in liquid withdrawal, coupled with the operation of adjacent injection wells, would have been expected to restore formation pressure – a restoration not observed in the field. This discrepancy is attributable to the dynamic lithological structure in the drainage area of well 3. In the geomechanically based simulation, the pressure drop led to a significant and irreversible reduction in permeability over large portions of the accumulation, thereby disrupting reservoir-to-well connectivity among production, injection, and aquifer regions. These results emphasize the importance of timely formation pressure maintenance in fractured reservoirs and highlight the critical role of accurately determining rock compressibility in flow simulations [32, 33].
Fig.5. Distribution of well 3 rock permeability across the section at the beginning of development (a) and after 10 years (b)
Fig.6. Dynamics of fluid flow rate and reservoir pressure of well 3
1 – fluid flow rate, m3/day; 2 – actual data; 3 – initial fluid flow model; 4 – fluid simulation model based on the geomechanical approach
Fig.7. Crossplot of the correspondence between the actual and estimated reservoir pressures in wells for the initial model (a) and the model modified using geomechanical modelling (b)
Figure 7 shows crossplots of actual versus estimated formation pressures for all wells. The geomechanically modified model shows a closer match to the actual pressures compared to the original model, without any deterioration in the prediction quality for well flow rates or water cuts.
Overall, our simulations indicate that the proposed concept effectively improves the fidelity of reservoir simulation models for hydrocarbon accumulations with fractured reservoirs.
Developing new approaches for reservoir simulation of fractured reservoirs is crucial for the design and development of hydrocarbon fields.
Although a discrete fracture network (DFN) model coupled with a dual-porosity/dual-permeability simulation remains the most comprehensive method for representing fracturing, it requires extensive parametric input.
Recently, three-dimensional geomechanical models have become increasingly popular because they can simulate the spatial distribution of principal stresses based on initial physicomechanical properties using the finite element method. While geomechanical models are often developed to support well drilling and hydraulic fracturing operations, their results can also be used to assess fracture aperture in fractured reservoirs and, by extension, to generate 3D permeability cubes.
Generalization of the fracturing study results enabled us to determine the general angle parameters of the fracture systems for the pay under consideration. Machine learning classification using angle stack cubes from 3D seismic surveys yielded a fracturing probability cube. A strong correlation was found between reservoir permeability and the rock state as determined by geomechanical modeling, enabling us to derive the full permeability tensor for the TFm pay.
The tNavigator software package implements a method for considering the effect of permeability reduction in flow simulations. The overall decrease in reservoir permeability is modelled using the built-in keywords of the simulator.
Reservoir simulation modelling using the described approaches has made it possible to significantly improve the convergence of formation pressure dynamics with an acceptable quality of convergence of production dynamics and water cut. Increasing the reliability of formation pressure forecasting allows planning inflow stimulation activities with a more predictable result, creating more realistic fracture properties when developing hydraulic fracturing designs, and optimizing the formation pressure maintenance system.
