Modeling reservoir permeability through estimating natural fracture distribution and properties

Permeability in earth models has three components: natural fracture, distinctive matrix permeability and porosity correlated matrix permeability. While matrix permeability is usually predictable from porosity and diagenesis effects, fracture permeability can be a highly uncertain parameter. The earth model permeability components are calibrated by determining fracture distribution and estimating fracture properties through an iterative optimization process. The calibration proceeds iteratively until the current estimate of calculated reservoir flow capacity is within acceptable accuracy limits to reservoir flow capacity indicated from production logging tool measurements.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to developing or forming models of subsurface reservoirs and estimating natural fracture distribution and properties to calibrate permeability for geological modeling of rock layers in subsurface reservoirs.

2. Description of the Related Art

Natural fractures present in subsurface formations are discontinuities representing a surface or zone of mechanical failure in the formation. Natural fractures have been formed over geological time as a result of movements and deformations within the subsurface rock over time. Natural fractures continue to be formed as a result of microseismic events which are slight tremors or movements in the earth's crust arising from various natural sources. Natural fractures are thus different in origin and nature from fractures induced in earth formations from the practice of hydraulic fracturing or fracking.

Natural fracture prediction is one of the more challenging problems in reservoir characterization. Fracture distributions are related to various factors such as intrinsic rock mechanics properties, as well as movements and deformation of the formation rock layers due to different tectonic stages to which the lithological formations are subjected through geological time. However, those parameters are usually unknown for the purposes of hydrocarbon exploration and development.

In highly complex geological environments due to different tectonics, understanding the natural fracture distributions and their properties can be important for conducting hydrocarbon exploration and development programs. A natural fracture can define a geological trap indicating possible flow passage for fluids in the reservoir. Natural fractures reduce the risk of unsuccessful results in drilling operations, and have an impact on reservoir management.

Hydrocarbons accumulate over geological time in a reservoir in the primary porous medium of the formation rock, and also in secondary porous media formed by natural fractures and in other areas of porosity such as vugs, caverns, and the like in the formation rock. Natural fractures as secondary storage mechanisms for hydrocarbon accumulation play an important role in some tight reservoir fields, enhancing the capacity to produce hydrocarbons from such reservoirs. Natural fractures enhance the permeability and connectivity between the primary porous media of the formation rock, and also support the flow of hydrocarbons into the wellbore. Natural fractures can also connect the porous and non-porous media of different rock layers of a reservoir in lower permeability conditions or situations.

Several techniques to build a natural fracture model have been described in the past, using different approaches. Examples include using seismic attributes from seismic exploration surveys for detection and distribution of natural fractures. Examples of this methodology are described in A. M. Zellou, T. Royer, G. C. Robinson, P. Zahuczki, A. Kirali, “Fractured Reservoir Characterization Using Post-Stack Seismic Attributes-Application to a Hungarian Reservoir”, Proc. 68th EAGE Conference and Exhibition, 2006; where post-stack seismic attributes were used as a basis for natural fracture prediction.

Other approaches for natural fracture prediction have been used for formation rock near a wellbore. These approaches have included fracture characterization based on measurements detected along the wellbore through borehole imaging (resistivity or sonic), measurements based on rock properties (porosity, density, etc.) obtained from rock core samples and well logs, and by fracture characterization along the wells. For areas of a reservoir distant from wellbores, this approach was often of limited applicability, since the measurement detection capabilities of well logging tools and the ability to obtain rock core samples has been limited to regions in the immediate vicinity of the wellbore.

Another conventional approach has been to design a reservoir fracture network in a deterministic way into a 3D geological model using conceptual models and validating the results with engineering data. This approach is based on manual mapping of natural fractures based on the numbers of wells intersecting fractures and seismic discontinuity attributes. The number of fractures could be increased as needed during a history match process to match the well test data. An example of this approach can be found in Rogers. S., Enachescu C., Trice R., and Buer K., 2007, “Integrating discrete fracture network models and pressure transient data for testing conceptual fracture models of the Valhall chalk reservoir”, Norway North Sea; Geological Society, London, Special Publication 2007, v. 270; p. 193-204. Other examples are in U.S. Pat. No. 7,565,278 (Li et al.); U.S. Pat. No. 9,390,204 (Bowen et al.); and U. S. Published Patent Application No. 2015/0276979 (Paradigm Sciences).

So far as is known, deterministic design of a reservoir fracture network relied on assumptions taking the fracture description from the wells and extrapolating those fracture descriptions into a grid model following a conceptual model. Such an approach was based on the reservoir being assumed to be one capable of representation as a simplistic model. Some assumptions for deterministic reservoir fracture network design were related to the natural fracture intensity (number of fracture/ft.). This deterministic method did not take into account fracture intensity indicated by borehole image logs or core descriptions from laboratory testing of formation core samples. Deterministic design assumptions also reduced the ability to predict or quantify natural fractures far from well locations. This additionally reduced reliability of reservoir fracture network model formed by deterministic design.

Fracture characterization and modeling of a reservoir has, so far as is known, been limited based on the acquired static data. The static data incorporates the geological features of the reservoir, including structure, stratigraphy, lithology and petrophysics of the subsurface formation rock.

Once a discrete fracture network has been identified, the accuracy of the fracture network model must be determined and confirmed. This has been done based on dynamic data, which are time dependent measures of flow and pressures in the reservoir wells under test conditions.

Previous techniques estimated permeability of the geocellular earth model considering a matrix component and fracture component. The fracture component was empirically calculated from theoretical correlations based on outcrop measurements relating fracture aperture to permeability. Multipliers were then identified to satisfy the pressure transient response or the production logging test (or PLT) flow profiles. The dynamic data has been used as input conditioning to alter data from which the discrete fracture network components have been determined. However, discrete fracture network fracture distribution and fracture permeability are independent physical conditions of a formation rock matrix.

If a well test showed high permeability in a region of interest, previous methods of fracture characterization and modeling indicated or assigned an area of increased permeability. However, so far as is known, previous methods did not differentiate whether the indicated high permeability was caused by fractures or high permeability streaks. With a model resulting from explicitly assigning a measured well test permeability around the zone of interest, it was not possible to determine whether the permeability resulted from a fracture, from vugs or from a high permeability matrix. Thus, the areal extent of the permeability phenomena could not be defined. This resulted in poor model predictability.

SUMMARY OF THE INVENTION

Briefly, the present invention provides a new and improved method of drilling a well in a subsurface geological structure to a location in a subsurface hydrocarbon reservoir indicated by a natural fracture network model of the reservoir. Reservoir parameters representing properties of the subsurface reservoir are obtained for processing in a data processing system. The natural fracture network model is then formed by processing the obtained reservoir parameters in the data processing system to identify fracture properties comprising the character, location, and stress conditions of natural fractures at locations in the subsurface hydrocarbon.

The processing determines fracture distribution based on the obtained reservoir parameters. A measure of estimated fracture properties of the subsurface reservoir is formed based on the determined fracture distribution. An estimated reservoir flow capacity is then obtained based on the formed estimate of reservoir fracture properties. An indicated flow capacity of the reservoir is also obtained. The estimated reservoir flow capacity is compared with indicated flow capacity of the reservoir.

If the estimated reservoir flow capacity is within acceptable limits of accuracy of the indicated flow capacity, a well is drilled in the subsurface geological structure to a location in the subsurface hydrocarbon reservoir based on the identified presence and extent of natural fractures in the subsurface geological structure. If not, an adjusted fracture distribution is determined and the formed estimate of reservoir fracture properties adjusted based on the adjusted fracture distribution. The steps of forming a measure of estimated fracture properties, obtaining an estimated reservoir flow capacity, and comparing the estimated reservoir flow capacity are then repeated.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

In the drawings,FIG.1is an isometric view in schematic form of subsurface geological structure S or formations in the earth at a location where a subsurface hydrocarbon reservoir R in the form of a hydrocarbon producing formation rock layer10is present. As shown inFIG.1, the hydrocarbon producing formation rock layer10is present beneath several other formation rock layers, such as indicates at12,14and16below the earth surface20. As indicated at22,24and26, exploratory or production wells have been drilled to penetrate the earth formations through wellbores as indicated at28and30.

FIG.2is a much enlarged view showing a schematic three-dimensional very small segment or portion32of the subsurface hydrocarbon producing formation rock layer10ofFIG.1. The segment32is formed as indicated at34of the primary porous medium of the formation rock. An irregular system of microscopic fractures36and small cavities or vugs38are typically present in the primary porous rock medium34. Natural fractures in reservoirs can also be present across a wide range of scale, ranging from those shown inFIG.2in the form of microfractures to extensive fractures or faults of thousands of meters. The vertical extent of natural fractures is often controlled by thin layers in the form of shale beds or laminations, or by weak layers of rock in carbonate sequences in the earth.

As shown schematically inFIG.3, a schematic diagram of the methodology of the present invention for reservoir hydrocarbon exploration, and in particularly the location and completion of wells for hydrocarbon production is illustrated schematically at100. As indicated at102inFIG.3, reservoir parameters and properties from a plurality of disciplines of earth science are obtained, assembled and stored in a data processing system D (FIG.7). As shown at102, the reservoir parameters include seismic attributes from seismic surveys as indicated at104; rock and mechanical properties from geological modeling as indicated at106; measures from structural restoration models as indicated at108; rock geological characterizations as indicated at110obtained from formation core samples and well logs performed in the wellbores such as32and34; and reservoir engineering measures obtained as indicated at112from production measures and reservoir simulations of the reservoir layer10.

As indicated at W, and as shown in more detail inFIGS.4through10, the assembled reservoir parameters from step102are subjected to reservoir modeling and with natural fracture distribution and properties determination. The drilling of additional wellbores of the types already present, as indicated schematically at32and34, is then performed as shown120. The drilling during step120is at locations indicated appropriate by the models resulting from processing steps according to the workflow W. Drilling during step120is thus directed to regions of the reservoir layer10where fractures of the types conducive to increased production are likely to be present. Drilling during step120is also enhanced by drilling to avoid formations or layers regions not indicated to be hydraulically conductive. Drilling during step120is also improved by avoiding areas indicated to contain fractures likely to cause complications in drilling operations or otherwise adversely impact drilling operations. With the present invention, wells are drilled based on the fracture modeling of natural fracture characteristics of subsurface reservoir formations. The fracture models reflect the influence of several measures of formation parameters obtained from several disciplines.

The processing according to workflow W is shown schematically inFIG.4designates generally a workflow for reservoir modeling of natural fracture distribution and properties according to the present invention. Permeability in earth models has three components: natural fractures, distinctive matrix permeability, and porosity correlated matrix permeability. The distinctive matrix permeability is normally predictable and can be obtained from porosity and diagenesis effects. However, determination of fracture permeability to a useful degree of confidence has, so far as is known, been difficult.

The present invention provides measures of determined fracture distribution for a geocellular model of a subsurface reservoir, and provides estimated fracture properties of the model. This is done, as will be described, by determining fracture distribution and estimating fracture properties through an iterative optimization process. This process progresses iteratively until reservoir flow capacity of the model, based on estimated fracture properties, differs from actual measured reservoir flow capacity by less than a specified accuracy or convergence limit. The present invention provides a discrete fracture network scaled to a geocellular earth model of the subsurface reservoir, and wells matrix calibrating coefficients for use in reservoir simulation and in estimating reservoir production.

Methodology

The present invention provides a methodology to validate a discrete fracture network in the reservoir, and forms estimates of earth model permeability based on the parametrization of multiple stochastic discrete fracture network (DFN) realizations or postulated fracture distributions, including whether the fractures are hydraulically conductive of fluids.

A hierarchy is established to calculate the equivalent permeability for the three components, where fractures have the major impact for the fluid flow movement, followed by high permeability streaks or HPS, then conventional matrix porosity. Using this hierarchy, the flow capacity is calculated for each component and optimized using the reservoir dynamic response to production logging tests.

As will be explained, the determination of fracture permeability is based on stress conditions on the formation rock. The present invention then determines through optimization processing an optimum set of permeability component parameters to satisfy a measured dynamic reservoir pressure and flow response to a production logging tool test.

Petrophysical Properties and Principles

The natural fractures prediction according to the present invention forms a geomechanical model of the subsurface rock. Formation of the geomechanical model includes in the processing an in situ stress regime and the effects of paleo-stress deformation accumulated over geological time. In the context of the present invention, the in situ stress regime is a condition where the stress field is unperturbed or is in equilibrium without any production or influences of perforated wells. The natural fracture system is closely related with current and past deformation due to the stress variation through the geological time. Different types of fractures can be formed in the formation rock during those deformation episodes.

As shown schematically inFIG.7, the methodology of the workflow W is organized into two components or stages: a Data Input and Conditioning Stage or Component C and a Process Engine Stage or Component E.

Data and Input Conditioning Stage C

The Data and Input Conditioning Stage C receives input measurements from laboratory tests, well logs, well measurements and performance data for present invention. The measurements are of three types or categories, which are: reservoir dynamic measurements, as indicated by a Reservoir Dynamic Measurements Module220(FIG.5); matrix and high permeability streaks (HPS) as indicated by a Matrix and HPS Module230(FIG.6); and geomechanical drivers to constrain discrete fracture network (DFN) realizations, as indicated by a Geomechanical Fracture Drivers Module240(FIG.7).

Reservoir Dynamic Measurements Module220

The reservoir dynamic data are acquired by field measurements in the wells, these measurements quantified the total or segregate fluid production, total flow capacity, damage factor, etc. In most of the case this kind of test are unable to evaluate the natural fracture response and distinguish it from matrix or HPS features. However, the total flow capacity (KH) can be estimated as well as the contribution of individual intervals. These measurements are obtained from two kind of test the PLT (Production Log Tool)222and PTA (Pressure Transient Analysis)224.

Production Log Tool (PLT)222

Production logging tools in well boreholes in a reservoir measure as indicated at222the nature and behavior of fluids and fluid distribution in or around the borehole during production or injection. Production logs are used to analyze dynamic well performance and productivity or injectivity, in particular the production or injection flow profile as function of depth in a well bore. The production logging tool measurement data222is provided as indicated at223as an input to a dynamic reservoir model permeability optimization module260(FIG.9). The production logging tool measurement data222is provided as indicated at223as an input to a dynamic reservoir model permeability optimization module260(FIGS.4and9).

Pressure Transient Analysis (PTA)224

The pressure transient analysis performed as indicated at224determines a measure of indicated or estimated actual flow capacity from producing formations in the reservoir. The analysis is based on measurements of the reservoir pressure changes over time during what are known as pressure transient tests. In such testing, a limited amount of fluid flows from the formation being tested, and pressure at the formation is monitored over time. Then, the well is closed and the pressure monitored while the fluid within the formation equilibrates. Analysis of these pressure changes provides information on the size and shape of the reservoir, as well as its ability to produce fluids. The flow total flow capacity (KHPrA) of a well is a main input variable utilized according to the present invention workflow. The indicated or estimated actual flow capacity results from pressure transient analysis are provided as indicated at225to the dynamic reservoir model permeability optimization module260. The indicated actual flow capacity results from pressure test analysis224are provided as shown at225to the dynamic reservoir model permeability optimization module260.

Matrix and HPS Permeability230

Matrix permeability232(KFacim=KMatrix), can be calculated and modeled from petrophysical workflows using a combined core test and wireline log data. The matrix permeability must be modeled into a conventional 3D geocellular grid, using algorithms for correlations and distribution.

There are typically high permeability streaks (HPS) in a producing formation matrix. The high permeability streaks are separate from natural fractures and are features related mainly to carbonate environments. High permeability streaks correspond to intervals with high permeability due to the presence of vuggy porosity or other diagenetic features in the rock matrix. The presence of HPS intervals can be detected partially by the conventional permeability workflow for matrix using a comprehensive characterization of the rock types. However, in some cases those intervals cannot be quantified correctly due to bad borehole conditions and may thus be underestimated. In such cases, additional corrections may be provided, such as those from what is known as a MHPs parameter. The matrix permeability model formed at232is provided as shown at233to the dynamic reservoir model permeability optimization module260.

Petrophysical rock Type (PRT)234inputs are data relating to the petrophysical rock types. The PRT input data234are computed based on core-plugs test laboratory measurements. The relationships between the permeability and porosity for the rock matrix are subdivided by using pore throat size distribution and a known theoretical model such as the Winland (R35) method. The output of the petrophysical rock type process234is data indicating the quality of the rock in terms of its permeability in order to quantify the High Permeability Streaks (HPS) for the 3D geocellular model.

Porosity236: A matrix porosity model236is computed based on results obtained from a porosity well log derived from petrophysical interpretation. The porosity well log is scaled up to the 3D grid resolution. Geostatistical algorithms and correlations are applied to determine porosity of formation rock at locations between the wells and as a result the porosity is modeled. The matrix porosity model formed as indicated at236is provided to the dynamic reservoir model permeability optimization model260as indicated at237.

Discrete fracture network realizations formed according to the present invention are constrained by geomechanical drivers. The parametrization of the main variable to constrain the fracture presence or position into a 3D geo-cellular grid includes several physical characteristics of the rock matrix, typically including: fracture density, fracture length, fracture orientation and geometry. The physical characteristics are obtained according to the present invention based on input obtained regarding matrix and provided the data processing system D in the form of controlling by Borehole Image Interpretation (BHI), brittleness property248, paleo-stress analysis and an “in-situ” stress regime244. The fracture aperture and permeability are determined by critical stress analysis, which relates the stress distribution and the fracture planes.

Critical Stress Analysis242

During step242, a fracture model is formed based on critical stress analysis. Conventional methods to calculate fracture aperture (such as outcrops observations, core studies, and borehole image logs) have, so far as is known, generally produced heterogeneous apertures in the rock. Other methods have assumed constant aperture presence for a fracture model.

FIG.11Ais a schematic diagram of a fracture model F-1of an in situ heterogeneous stress field with a constant aperture Pfin a representative segment of a subsurface formation. However, with the present invention it has been found that a conventional model such asFIG.11Ais not directly representative of stress conditions during actual flow through a fracture network. Analysis of subsurface data indicates not all fractures contribute to flow. According to the present invention, a formation fracture model is formed, as will be described, at242. The formation fracture model so formed has only a determined portion of the discrete fraction network hydraulically open for passage of flow. Determination of such a portion and its distribution in the rock is based on whether the fractures are critically stressed. In determining critical stress, a physical phenomenon known as the Coulomb friction criterion is applicable.

The Coulomb criterion depends on the stress magnitude and the orientation of the fracture in the in situ heterogeneous stress field present in a formation. The fracture orientation with respect to the stress directions has significant impact on determination of normal and shear stresses on a fracture plane. When shear stress exceeds shear stiffness, shearing causes dilation keeps the fracture hydraulically open. Fractures in this stress state are referred to as reactivated or critically stressed. Fractures in this stress state are according to the present invention defined as critically stressed fractures in calculating the fracture aperture as a function of the shearing dilatations, as will be described.

Critical stress analysis is a function of normal stress σn, shear stress τ and fluid pressure. In the example fracture model F-1shown inFIG.11A, normal stress σnis expressed according to Equation (1) as follows:
σn=0.5*(σ1+σ3)+0.5*(σ1+σ3)*Cos 2θ  (Equation 1)

The shear stress τ is expressed according to Equation (2) as follows:
τ=0.5*(σ1−σ3)*Sin 2θ  (Equation 2)

In Equations (1) and (2), σ1and σ3are the maximum and minimum horizontal stresses in a horizontal plane-strain cross section of the fracture network Pf. θ is the angle between the plane of normal stress on and the direction of maximum stress σ1as shown inFIG.11A. Further discussions are contained inReservoir Geomechanics, Afark D. Zoback, Cambridge iUniversity Press, UK, 2007.

The relation between normal stress σn, shear stress τ and coefficient of friction ϕ in a rock matrix is represented graphically in what are known as Mohr diagrams or circles. An example Mohr diagram M (FIG.11B) indicates rock coefficient of friction ϕ as a function of shear stress τ and effective normal stress σnfor the heterogeneous stress conditions indicated schematically inFIG.11A.

In the context of the present invention, the Mohr's diagrams are graphical representations in two or three dimensions of stress conditions in a rock mass at different planes oriented as functions of orientation angle θiof planes passing through a point of interest in the rock. The Mohr's circles permit determination at the point of interest of principal normal stresses σmaxand σmin, maximum and minimum shear stresses τmax, τminas well as the orientation of the principal planes. Details of Mohr's circles are explained, for example, in Mohr's Circle. Mohr's Circle, Wikipedia. April, 2019, https://en.wikipedia.org.Mohr%27s_circle. As will be explained, the Mohr's circles inFIG.14Bof a critical stress fracture model F-2inFIG.14Ainclude a plot of rock coefficient of friction ϕ to indicate the physical phenomenon of critical stress.

Critical stress concept criteria are used during critical stress analysis processing242in the methodology of the present invention to obtain natural fracture properties. This is done governed by the Coulomb criterion, which depends on the stress magnitude and the orientation of the fracture plane with respect to the “In situ” stress orientation. The orientation impacts the normal and shear stresses on the fracture plane. The results of the critical analysis states242are provided as indicated at243as an input to stochastic discrete fracture network model250.

When shear stress exceeds shear stiffness, shearing causes dilation and keeps the fracture hydraulically open, as described in Rogers, S. F. (2003),Critical stress-related permeability in fractured rocks, Geol. Soc. London, Spec. Publ.,209(1), 7-16, doi:10.144/GSL. SP.2003.209.01.02. Fractures in this example critical stress state of fracture model F-2are shown at340(FIG.14A) are referred to be reactivated or critically stressed in contrast to non-critical stress fractures342. This geomechanical feature is described in Barton, C. A., et al. (1995),Fluid flow along potentially active faults in crystalline rock, Geology,23(8), 683-686, doi:10.1130/0091-7613(1995)023<0683:FFAPAF>2.3. CO;2, and in Rogers, S. F. (2003), previously cited.

The “in-situ” stress regime of natural fracture distribution is preferably modeled using FEM software (Finite Element Model) techniques, which predicts a stress/strain tensor regime using mechanical boundary elements. FEM methods use geomechanics simulations to converge a proper solution under certain boundary stress conditions. Maximum principal horizontal stress model and magnitude are obtained by this methodology for each cell into the 3D Grid geo-cellular model as described in Herwanger, J., et al.: “Seismic Geomechanics—How to Build and Calibrate Geomechanical Models using3D and4D Seismic Data”,1Edn., EAGE Publications b.v., Houten, 181 pp., 2011. There are several available software applications to model stress regimes, such as Visage finite element geomechanics simulator from Schlumberger, Ltd., or Abaqus FEA finite element analysis suite from ABAQUS, Inc.

The Borehole image information obtained from borehole image logging tools during step110(FIG.3) is interpreted during the Geomechanical Fracture Driver processing246to provide natural fracture descriptions, including natural fracture type, dip angle, dip azimuth and intensity at well level. For natural fractures, the apparent aperture can also be estimated by using a normalized image through deep resistivity response. An apparent fracture aperture can be calibrated using real core-plug measurements at reservoir conditions. The borehole image interpretation data so formed during this stage is provided as input data as indicated at247to the Stochastic Discrete Fracture Network250for further processing. The fracture description is used to model the paleo-stress and geomechanical restoration.

The geomechanical restoration process is used to calculate stress and strain paleo-stress deformation by analyzing each geological tectonic episode. This analysis should preferably distinguish between the fractures created by folding process or by faulting process, creating possible strain/stress deformation for each process. Paleo-stress analysis regarding the fractures folding relation may be modeled using geomechanical restoration functionality such as that available in geological structural modeling analytical methodologies such that available as Kine3D from Emerson Paradigm Holding LLC., or that provided as Move Suite from Petroleum Experts, Inc. The faulting response may be modeled using boundary element methods (BEM) which are incorporate into available processing methods, such as Petrel software iBEM3D as described in “Adaptive cross-approximation applied to the solution of system of equations and post-processing for 3D elastostatic problems using the boundary element method,” Engineering Analysis with Boundary Elements 34 (2010), p. 483-491, F. Maerten.

A rock brittleness property is differentiated during step248. When subjected to stress, brittle rock breaks without significant plastic deformation. Brittle materials absorb relatively little energy prior to fracture, even those of high strength. In a complex heterogeneous fracture rock mass, the brittleness property is modeled using neuronal network classification, taking as inputs the elastic properties and stress regime producing mechanical facies. Those mechanical facies preferably should have some proportional relation with the distribution of natural fractures. This correlation can be evaluated using histograms filtered by density fracture. The rock brittleness property data from states248is provided as shown at249as an input to the stochastic discrete fracture network250.

Process Engine Stage E

The discrete fracture network or DFN processing250is a module which according to the present invention forms a computational stochastic discrete fracture network model on a stochastic or statistically random basis from inputs provided. The stochastic discrete fracture network250receives as inputs the borehole image information from geo-mechanics fracture driver246and rock brittleness property data from the module248. Operation of the stochastic discrete fracture network250begins at a start step252to begin formation of a new realization of a computerized discrete fracture network model. Next, a discrete fracture network is formed representing parameters of fracture length and intensity during step254based on the received inputs. Step254may be performed, for example, according to the methodology of commonly owned U.S. patent application Ser. No. 15/704,236, filed Sep. 14, 2017, “Subsurface Reservoir Model with 3D Natural Fractures Prediction.” Step256follows during which new values for a fracture aperture are formed in terms of area, porosity ϕ, shear stress τ and effective normal stress σn. Step256is based on the results of the critical stress analysis242. The stochastic discrete fracture network model so formed during step256represents geometrical properties of each individual fracture (e.g. orientation, size, position, shape and aperture), and the topological relationships between individual fractures and fracture sets.

A number of discrete fracture network realizations or postulated fracture distributions are formed during processing according to the present invention controlled by several sets of fracture parameters. The parameters for this purpose include fracture length, density, orientation, geometry, aperture and permeability. The individual discrete fracture network realizations so formed are distributed into a 3D geocellular grid to produce reliable stochastic fracture realizations, existing platform such as the Petrel™ fracture modeling package or other suitable commercially available software adapted to populate natural fractures.

The initial four parameters are distributed based on mechanical earth model properties such as borehole image interpretation, mechanical layering, paleo-stress inversion and the in-situ stress regime. The fracture aperture distribution is controlled according to the present invention by whether a fracture is at critical or non-critical state. This in turn controls whether the fracture to be hydraulically conductive or not, as has been described.

Stochastic fracture network realizations are formed based on continuous property model, such as fracture density, orientations and geometry. The geomechanics fracture driver modeled in the steps of: borehole image interpretation110, rock brittleness property248and stress regime model242are provided as inputs to the stochastic discrete fracture network250for production of a probabilistic fracture realization for the model being formed.

The stress regime predicted for the 3D grid model during previous step244is used in processing by stochastic discrete fracture network250to apply the Coulomb failure criteria to the fracture planes generated in the DFN models. As will be explained, application of Coulomb failure criteria results in the hydraulic closing of all fractures in the rock matrix which are not within the critically stressed orientation domain.

The fracture aperture Pfis based on normal closure and shearing dilatation. Only near-critically-oriented fractures which are oriented to be subject to stress distinctions which are near critical can dilate. Shear dilation occurs only partially in the natural fractures present in the rock matrix. Other fractures in the rock matrix which are not subject to near critical stress do no exhibit dilation.

The physical phenomenon of near critical stress and its determination to model stress-dependent permeability for fractures are discussed in “Stress-Dependent Permeability of Fractured Rock Masses: A Numerical Study”, Ki-Bok Min et al.,International Journal of Rock Mechanics and Mining Sciences, Volume 41, Issue 7, 2004, Pages 1191-1210, ISSN 1365-1609. The present invention further incorporates the effects of both normal closure and shear dilation of fractures as described for example in “A continuum model for coupled stress and fluid flow in discrete fracture networks” (Quan Gan et al.,Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:43-61; and “The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks” (Qinghua Lei et al.,Computers and Geotechnics journal,2016.)

FIGS.12and13illustrate basic mechanisms of aperture changes resulting from two basic fracture deformations as observed in numerous laboratory experiments and the states of stress on fractures. The change in fracture aperture occurs from two basic mechanisms.FIG.12illustrates changes in fracture aperture as a result of normal stress-induced closure or opening.FIG.13illustrates changes in fracture aperture as a result of shear stress-induced dilation.

As shown at332inFIG.12, in normal stress-displacement relation the response shows well-known hyperbolic behavior of stiffer fractures with increasing compressive normal stresses. The normal deformation of a fracture is more sensitive at lower magnitudes of normal stresses.

Shear dilation occurs as a result of overriding asperities of two rough fracture surfaces and may reach a stationary value with increasing fracture shearing, as shown at334inFIG.13. Such shear dilation of fractures is important, since its magnitude can be significantly larger than the normal stress-induced closure or opening. The model formed according to the present invention considers permeability changes from normal and shear stresses as being independent of each other.

FIG.14Ais a schematic diagram of a fracture model F-2and14B is a Mohr's diagram provided as examples illustrative of the interrelation of fracture aperture orientation angle Θito shear stress τθi; normal stress σnθi; and rock coefficient of friction ϕ. As indicated schematically inFIG.14Bfor a fracture orientation angle θithe shear stress τθiand normal stress σnθiconditions define whether the rock is or is not under critical stress according to the rock coefficient of friction ϕ.

As indicated inFIG.14A, fractures340are in critical stress. The stress conditions σnθ1and τθ1are indicated for fracture orientation angle Θ1in the Mohr's diagram ofFIG.14B. The critical stress conditions indicated in Mohr's diagram M-2(FIG.14B) for the fractures340are greater than the coefficient of friction (ϕ). This indicates that under the formation fracture angle and critical stress condition shown inFIG.14A, the formation rock in the fracture model F-2is hydraulically conductive.

The equivalent aperture dilatation of normal closure and shear dilatation phenomena can be expressed according to on empirical Equations (3) and (4):

In Equations (3) and (4), τ is the shear stress and σnis the gradient of normal stress acting in a single plane; φ is the friction angle; α is a coefficient related to the fracture closure; c is a stress coefficient for dilatation; b0is the residual aperture; dmaxis maxima dilatation and bmaxis the maximum mechanical deformable mechanical aperture.

Using the previous aperture distribution, a standard cubic law function is used as expression for the stress-dependent permeability (Equation 5). The stress-dependent permeability so expressed incorporates the effects of both normal closure and shear dilation of fractures through the aperture distribution based on critical stress:
K≈B*(Aperture)2(Equation 5)

The B parameter for Equation (5) is calibrated using multiple realizations and comparing the resultant realizations with the dynamic well test response shown inFIGS.8and9. The B parameter is a constant value parameter for the intrinsic permeability multiplier and ranges are initially assumed from theoretical correlations. The B parameter is then optimized through the workflow.

Fluid Pathways

Fracture networks usually serve as the major pathways for fluid transport in subsurface rocks, especially if the matrix is almost impermeable compared to the fractures. The partitioning of fluid flow within a fracture population relies on the spatial connectivity of fracture geometries and the transmissivity of individual fractures, both of which can be affected by the geomechanical conditions.

FIG.15Ais a schematic diagram of a formation model F-3of fluid pathways in formation rock in a segment of a subsurface formation.FIG.15Bis a Mohr's diagram indicating critical stress conditions in certain of the fluid pathways in the formation model F-3. The stress conditions are indicated for the fluid pathways in fracture model in the Mohr's diagram ofFIG.15B. The critical stress conditions indicated in Mohr's diagram M-3(FIG.15B) for the certain fluid pathways are greater than the coefficient of friction (ϕ). This indicates that under the fluid pathway angle and critical stress condition shown inFIG.15A, the formation rock in the formation model F-3is hydraulically conductive.

The flow paths350based on critical stress (FIG.15A) which are hydraulically conductive are indicated in red, while flow paths352indicated in blue do not conduct fluids. The presence or absence of critical stress which affects hydraulic conductivity in natural fractures is a relationship as follows according to Equation (6) from Zoback,Reservoir Geomechanics, previously cited:
Critical fractures=(τ−σn*Tan(φ))≥0  (6)

In the dynamic permeability optimization processing module260, as indicated at step262static reservoir data in the form of flow capacity of the reservoir matrix Kfacimstored in memory during step32, and rock type (PRT) stored in memory during step34are received and stored for further processing.

As indicated at264, a sequence of discrete fracture realizations are generated by the discrete fracture network realization250are received and stored as indicated at step266. The discrete fracture network is also during step266upscaled from the 3-dimensional object planes space to the geocellular gridded model. The natural fractures are represented into a 3D grid model as planes in the explicit way which contains specific aperture and intrinsic permeability, in order to convert the fracture planes in a Geocellular properties upscaling process is made using a software package such as Petrel™, Fracflow, or other suitable upscaling methodology.

The upscaling is performed to form fracture porosity distribution, transfer functions coefficients, and fracture permeability measures which are stored for further processing. The stored tensor effective permeability (Ki, Kj, Kk) parameter measures for the current discrete fracture network realization are provided as indicated at step268for determination of the flow capacity KHFract-simulateof fractures in the rock matrix.

The stored reservoir data Kfacimfrom step262are provided as indicated at step270for comparison with the flow capacity dynamic data from production logging tool well tests. If the comparison during step270indicates that the well flow capacity KHPTAis less than the flow capacity KHFract-simulateof the rock matrix, an indicator is provided as step272that a problem exists in the data being provided for processing.

If the comparison during step270indicates no problem is present in the data, processing proceeds to step274. During step274, the determined flow capacity KHFract-simulatefor the fractures in the discrete fracture network realization in the sequence currently being processed is examined. This is done to detect null values in the flow capacity KHFract-simulateto identify fracture free wells as contrasted with wells which are intersecting fracture planes and exhibit a non-null determined flow capacity KHFract-simulate. If a null value is indicated during step274for the flow capacity KHFract-simulatefor a realization of a well, the identification of that well realization stored as indicated at276in a well list in memory of the data processing system D.

Wells indicated to be not intersecting fractures are treated during step278by conventional matrix calibration, a global multiplier is identified based on pressure transient analysis (PTA) results, and then propagated through individual grid block according to the production logging tool (PLT) flow profile distribution.

If an actual, non-null value of flow capacity is indicated in step274, processing proceeds to step280for classification processing relating to processing for wells which are intersecting fractures. During step280the presence of high permeability streaks is evaluated using the PRT corresponding to HPS code and if a particular well is indicated to contain a high permeability streak then this can be checked through the expression PRT=HPS. If presence of high permeability streaks is indicated for a well, the well is identified during step282in the database memory of the data processing system D as having a high permeability streak present.

If a well is indicated that during step280as having no high permeability streak, the well is so identified in memory of the data processing system D as indicated at284. Processing proceeds to step286to determine if the well is one which has had a production logging tool (PLT) well test performed and data obtained. For wells in which there is no PLT log data obtained, during step288, the difference E0 between the total measured flow capacity KHPTAand each of the predicted flow capacity KHFacimand the flow capacity for fracture (KHFract-simulate) is determined according to the relationship:
E0=KHPTA−KHFacim−KHFract-simulate(7)

When PLT log data is available from the well, the processing stage289follows during which parameters are defined for flow capacity determination during subsequent processing. The parameters and their designations are indicated in the drawings as follows:
HMS=1−HFS
PLTF=ΣPLTn
PLTM=ΣPLTMI
KHPTAF=PLTFKHPTA
KHPTAM=KHPTA−KHPTAF

The difference E1 between the fraction of the total flow capacity attributed to fractures KHPTAFand the accumulated flow capacity Σ KFSiHFSiis determined during step290through the relationship:
E1=(KHPTAF−ΣKFSiHFSi)  (8)

In Equation (8), the parameter KHPTAFrepresents the fraction of the total flow capacity attributed to fractures and ΣKFSiHFSirepresents the accumulated flow capacity for the PLT points being simulated (intersecting fractures).

During step292, for well depth intervals between a well depth or PLT point at which production logging tool data are obtained in order to determine the flow contribution at a particular depth in the well, the difference Ei between the portion of flow capacity assigned from the total flow capacity to each PLT point PLTFi−KHPTAFand the simulated flow capacity KFSiHFSiis represented by the following expression:
Ei=(KFSiHFSi−PLTFi·KHPTAF)  (9)
where the parameter PLTFirepresents the percentage of flow contribution at a single point of well depth, and KFSiHFSi represents the flow capacity contribution for the PLT points for fractures simulated; as shown in Figure ?.

Step294follows, during which the difference Ei for each data point determined during step292are assembled into an error matrix. Processing proceeds to a minimization processing procedure as indicated at step300(FIGS.9and19).

The minimization processing procedure is performed for each realization, and thus E0, E1, Ei are estimated. These are considered the error functions to be optimized (minimized), and minimization processing300is performed to select the best parameters to fit the dynamic of the formation rock to the results of the production logging tests.

Performance of minimization processing according to step300minimize the differences for each realization stored in the error matrix during step294as error factors for each of the wells across all their connecting grid blocks in the geocellular model.

The minimization processing300generates a population of realizations varying the parameters controlling the fracture distribution and properties as defined by the stochastic discrete fracture network processing module250. The number of the realizations in a single population is user defined.

The minimization processing workflow during step300(FIG.16) determines the error factors and populates the error matrix for each realization. During step302, fracture densities, orientation and geometrics from the geomechanics fracture driver246are received as inputs. Step304is performed to determine an error value between the measured flow capacity for the well and the flow capacity determined for each of the realizations. The determined error values are populated in an error matrix during step306and ranked during step308based on the different grid block fit factor in the error matrix. During step310a convergence test is performed to determine if the current estimate of calculated reservoir flow capacity for each realization is within user specified acceptable accuracy limits to actual measured reservoir flow capacity. If not, during step312a realization count is incremented and as indicated at313processing returns for formation of a new realization of a stochastic discrete fracture network by stochastic discrete fracture network processing250.

This process continues iteratively until reservoir flow capacity of the model based on estimated fracture properties is within an acceptable degree accuracy with actual measured reservoir flow capacity. The degree of accuracy is user specified and is based on considerations of computer processing time requirements, prevention of iterative looping, and reservoir engineering requirements.

If during convergence step310each realization is within user specified acceptable accuracy limits to actual measured reservoir flow capacity the determined fracture model is stored as indicated at314in memory of the data processing system and provided for output analysis in connection with drilling wells during step120(FIG.1).

Once the minimization processing is satisfactorily completed, and the final configuration for the fracture distribution so determined is stored in the data processing system D, an additional high permeability streak or HPS calibration indicated at320inFIG.10is performed. The high permeability streak calibration320is performed for those wells which have been indicated during step282(FIG.9) as passing through high permeability streaks. The difference between flow capacity from PTA, KHPTA, flow capacity equation for HPS is calculated during step322through the relationship:
KHHPS=KHPTA−Kfacim(H−HHPS)−KFract simulateHFS(10)

The relationship of Equation 10 determined during step322indicates the difference between the flow capacity KHPTA, the flow capacity from the fracture simulated high permeability streak HPS, and the flow capacity from the matrix. Step324is performed to determine whether KHHPSis less than zero then, conventional correction is prorated during step326through the flow profile. If the KHHPSis greater than zero, then a multiplier factor is formed during step328and applied to attribute this difference in flow capacity to the HHPSlayers.

Data Processing System D

As illustrated inFIG.17, the data processing system D includes a computer400having a master node processor402and memory404coupled to the processor402to store operating instructions, control information and database records therein. The data processing system D is preferably a multicore processor with nodes such as those from Intel Corporation or Advanced Micro Devices (AMD), or an HPC Linux cluster computer. The data processing system D may also be a mainframe computer of any conventional type of suitable processing capacity such as those available from International Business Machines (IBM) of Armonk, N.Y. or other source. The data processing system D may in cases also be a computer of any conventional type of suitable processing capacity, such as a personal computer, laptop computer, or any other suitable processing apparatus. It should thus be understood that a number of commercially available data processing systems and types of computers may be used for this purpose.

The computer400is accessible to operators or users through user interface406and are available for displaying output data or records of processing results obtained according to the present invention with an output graphic user display408. The output display408includes components such as a printer and an output display screen capable of providing printed output information or visible displays in the form of graphs, data sheets, graphical images, data plots and the like as output records or images.

The user interface406of computer400also includes a suitable user input device or input/output control unit410to provide a user access to control or access information and database records and operate the computer400. Data processing system D further includes a database of data stored in computer memory, which may be internal memory404, or an external, networked, or non-networked memory as indicated at416in an associated database418in a server420.

The data processing system D includes program code422stored in non-transitory memory404of the computer400. The program code422according to the present invention is in the form of computer operable instructions causing the data processor402to form subsurface reservoir 3D geocellular models including natural fracture presence, distribution, and flow parameter properties according to the present invention in the manner set forth.

It should be noted that program code422may be in the form of microcode, programs, routines, or symbolic computer operable languages capable of providing a specific set of ordered operations controlling the functioning of the data processing system D and direct its operation. The instructions of program code422may be stored in memory404of the data processing system D, or on computer diskette, magnetic tape, conventional hard disk drive, electronic read-only memory, optical storage device, or other appropriate data storage device having a computer usable non-transitory medium stored thereon. Program code422may also be contained on a data storage device such as server420as a non-transitory computer readable medium, as shown.

The data processing system D may be comprised of a single CPU, or a computer cluster as shown inFIG.17, including computer memory and other hardware to make it possible to manipulate data and obtain output data from input data. A cluster is a collection of computers, referred to as nodes, connected via a network. Usually a cluster has one or two head nodes or master nodes402used to synchronize the activities of the other nodes, referred to as processing nodes424. The processing nodes424each execute the same computer program and work independently on different segments of the grid which represents the reservoir.

The present invention through fracture distribution optimization process identifies individual measures of both a predictive fracture model and distinctive fracture properties from matrix properties. The optimization process does not use the dynamic data as a conditioning input but as match criteria. This provides a high quality product with predictive capacity of the unknown occurrence or properties of fracture, this is a superior quality than simply reproducing the actual measured data. This will provide an early indication of the type and nature of fractures to be encountered which in turn permit adjustments and modifications of the well path or well completion method being performed. The present invention is thus integrated into a practical application in that it solves a technological problem by allowing control of drilling operations to take into account conditions of favorable presence of increased permeability in the formations of interest.

The present invention differs from previously used methods in that dynamic data from reservoir production logging tests is not used as input parameters or measures to determine fracture occurrence or permeability of the formation rock. The dynamic data is instead used as a component of an error function in an optimization processing to determine the predicted fracture model and matrix properties.

The error function formed for optimization indicates which realization of predictive fracture model and distinctive fracture properties for the formation provides a reservoir flow capacity of the model based on estimated fracture properties within an acceptable degree accuracy with actual measured reservoir flow capacity.

The present invention solves thus solves technological problem. The present invention provides measures to overcome this problem through separately determining each of natural fracture occurrence and distribution in formations for a geocellular model; formation permeability in high permeability streaks; and distinctive super permeable matrix properties.

The invention has been sufficiently described so that a person with average knowledge in the matter may reproduce and obtain the results mentioned in the invention herein Nonetheless, any skilled person in the field of technique, subject of the invention herein, may carry out modifications not described in the request herein, to apply these modifications to a determined structure, or in the manufacturing process of the same, requires the claimed matter in the following claims; such structures shall be covered within the scope of the invention.

It should be noted and understood that there can be improvements and modifications made of the present invention described in detail above without departing from the spirit or scope of the invention as set forth in the accompanying claims.