WO2019238451A1 - A method and a system for modelling and simulating a fractured geological structure - Google Patents

A method and a system for modelling and simulating a fractured geological structure Download PDF

Info

Publication number
WO2019238451A1
WO2019238451A1 PCT/EP2019/064318 EP2019064318W WO2019238451A1 WO 2019238451 A1 WO2019238451 A1 WO 2019238451A1 EP 2019064318 W EP2019064318 W EP 2019064318W WO 2019238451 A1 WO2019238451 A1 WO 2019238451A1
Authority
WO
WIPO (PCT)
Prior art keywords
fracture
macrofracture
timestep
microfracture
macrofractures
Prior art date
Application number
PCT/EP2019/064318
Other languages
French (fr)
Inventor
Michael Welch
Mikael LÜTHJE
Original Assignee
Danmarks Tekniske Universitet
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Danmarks Tekniske Universitet filed Critical Danmarks Tekniske Universitet
Publication of WO2019238451A1 publication Critical patent/WO2019238451A1/en

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • E—FIXED CONSTRUCTIONS
    • E21—EARTH DRILLING; MINING
    • E21B—EARTH DRILLING, e.g. DEEP DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B43/00—Methods or apparatus for obtaining oil, gas, water, soluble or meltable materials or a slurry of minerals from wells
    • E21B43/25—Methods for stimulating production
    • E21B43/26—Methods for stimulating production by forming crevices or fractures
    • G01V20/00—
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/64—Geostructures, e.g. in 3D data cubes
    • G01V2210/646—Fractures

Definitions

  • the present invention relates to a method and a system for modelling and simulating a fractured geological structure using a dynamic geomechanical model comprising an implicit fracture model and an explicit discrete fracture network model.
  • geological and physical principles for the modelling and simulation of the geology are the same for all industries, although there may be some differences in the outputs desired depending on the decision to be taken or the action to be performed. At present, such geological models are especially prevalent in the oil and gas industry where they can be used to plan optimal hydrocarbon extraction strategies.
  • the geological structures of interest often include fractured reservoirs. If there are interbedded layers in the geological structure with contrasting mechanical properties, the fractures may often be confined within brittle layers sandwiched between more ductile layers, forming“layer-bound fractures” and the overall fracture population within each brittle layer can typically be subdivided into sets of approximately parallel fractures.
  • Such fracture networks will have a major impact on fluid flow through a fractured layer, especially when the rock mass has low permeability, as they will provide critical long-distance fluid conduits. They will also have a critical impact on the strength and elastic properties of the fractured layer as a whole (the“bulk rock” properties), potentially introducing planes of weakness and elastic anisotropy to the layer.
  • the conventional solution to this problem is to try to characterise the subsurface fracture network as best can be done from available data (typically borehole data) and then extrapolate this to build a stochastic fracture model.
  • stochastic models (“implicit” fracture models)
  • the fractures are not represented individually, but rather the fracture network is described statistically by one or more spatially variable properties, such as fracture density, fracture orientation, or fracture connectivity. These properties are then used to directly modify key bulk rock properties, such as porosity and permeability (for flow modelling) or failure strength and elastic moduli (for geomechanical modelling).
  • TM Transmissibility Multiplier
  • the properties that can be measured and extrapolated in this way are usually very limited.
  • the boreholes provide no data on the fracture sizes or connectivity, so it is not possible to determine the length distribution (or even the mean lengths) of the fractures, nor is it possible to determine how many fractures are interconnected and how many are isolated. Indeed, it is not even possible to determine the average number of fractures per unit volume (the“volumetric fracture density”, known as P30). These properties are all important in controlling the fracture permeability.
  • DFN Discrete Fracture Network
  • the aim is to generate a DFN, which has the same overall statistical properties as the real fracture network.
  • Key properties to be replicated are the fracture density, (both the volumetric fracture density (P30) and the mean linear fracture density (P32), the fracture size distribution, and the fracture connectivity.
  • a third approach is to generate a fracture model by simulating the growth of the fracture network dynamically, based on fundamental geomechanical principles.
  • the main advantage of this approach is that no a priori knowledge of the fracture density, size distribution, connectivity or any of the other statistical parameters that are so difficult to measure in the subsurface is required. Nor is it necessary to specify the fracture nucleation and propagation rates. All that is needed is a knowledge of the geomechanical properties and deformation history of the fractured layer. This information is usually readily available: the
  • geomechanical properties can be taken from mechanical tests on core samples, wireline log geophysical data or published data for similar lithologies, while the deformation history is often described in studies of the regional geology, or can be inferred from the geometry of larger-scale geological structures (e.g. major faults, folds or diapirs).
  • DFNs generated in this way will be more realistic than stochastic DFNs as they will by definition be consistent with the geomechanical processes of fracture development, and will hence honor the geology of the fractured layer and the larger-scale geological structure.
  • the statistical properties of the fracture network emerge from the simulation, unlike the stochastic modelling approaches where they are specified in advance. These simulations therefore generate new information, unlike the stochastic DFNs which simply reflecting the (usually unreliable) data used to create them.
  • this approach involves using some numerical geomechanical method to simulate the propagation and interaction of individual fractures and thereby generate an explicit DFN.
  • the finite element method is most commonly used, but other methods such as the boundary element method or the discrete element method have also been used for this purpose.
  • This approach has been used to model artificially induced fractures as well as natural fracture networks.
  • Numerical geomechanical methods allow detailed modelling of local stress anomalies that build up around fractures and faults, as well as local variation in mechanical properties (e.g. cemented nodules, stiffer or more ductile layers, or previous fractures). Thus, they can predict the effect that these local variations will have on the ultimate fracture geometry.
  • CCS Carbon Capture and Storage
  • the present invention relates to a method for modelling and simulating a fractured geological structure using a dynamic geomechanical model, which is discretised into geological layers, timesteps and gridblocks and comprises an implicit fracture model and an explicit discrete fracture network model, wherein the method comprises an implicit fracture modelling procedure for simulating the processes of fracture nucleation, propagation and interaction, based on fundamental geomechanical principles, to generate an implicit fracture model, in which model a fracture population is represented statistically by functions or values, and an explicit fracture modelling procedure for simulating the process of fracture propagation, interactions between parallel fractures and intersections of non-parallel fractures in an explicit discrete fracture network model, in which a fracture is represented as an individual object, based on data from the implicit fracture model relating to fracture nucleation and propagation (such as rates of fracture nucleation and growth), as well as incorporating the explicit discrete fracture network model into a geological simulator.
  • the method according to the present invention does not only give a more realistic representation of the fractures, honoring the geomechanics and the geology but also improves the accuracy of the resulting dynamic geomechanical model of the geological structure of interest and reduces the time taken to generate it.
  • the use of a combined implicit and explicit model makes it possible to simulate fracture networks through a full range of different fracture sizes scales and across large geological structures within a reasonable time and using commonly available input data.
  • the method according to the present invention is able to model both Mode 1 dilatant fractures and Mode 2 shear fractures, and to realistically recreate the geometry of the fracture network developed by the interactions between two sets of fractures, striking perpendicular to the minimum and maximum horizontal strain orientations. It is also capable of generating a purely implicit fracture model, in which the fracture sets are described by cumulative fracture density distribution functions, and can be used to generate large fracture models (e.g. across an entire geological structure) very quickly.
  • the method allows for the determination of the impact of the growing fracture network on the bulk rock elastic properties, and thus the changes in the in situ stress and the rate of fracture nucleation and propagation as the fracture network grows.
  • This dynamic data from the implicit model, combined with the rules governing fracture interaction, can be used to generate a geomechanically-based explicit DFN without simulation of the propagation of each fracture numerically.
  • Implicit fracture models can be generated independently for every gridblock in a static geomodel, using mechanical property and horizontal strain data specific to that gridblock, and these results can be combined to generate a single, integrated explicit DFN spanning the entire geomodel.
  • This is a key development as it allows for the generation of large, geomechanically consistent DFNs that reflect the variability in geometry, mechanical properties and in situ stress across a large geological structure, for example a DFN containing all layer-bound fractures across an entire salt diapir, or a near wellbore DFN containing the larger microfractures that develop in the damage zone around a nearby fault. This cannot be done using previous techniques known within the art.
  • Other key advantages of this method over previous methods are:
  • the implicit component of this method can calculate cumulative density distribution functions for the volumetric fracture density P 3 o(S), the mean linear fracture density P 32 (S), the fracture porosity ( P(S) and the total stress shadow volume IJJ(S), for each fracture set, over a full range of scales from mm-scale microfractures to large layer-bound fractures and even small faults. Note that these cumulative density distribution functions will not necessarily follow a power-law relationship, even if the original“seed” microfracture populations do, due to the effects of differential growth rates, fracture interaction and changes in the in situ stress as the fracture network grows.
  • the method can be calibrated by comparing the implicit cumulative fracture density distribution functions against the cumulative fracture density distribution functions calculated for the explicit DFN. A close match will demonstrate that the implicit component of the method is correctly predicting the fracture populations, based on the premises concerning fracture nucleation and interaction.
  • the dynamic geomechanical model comprises one or more microfractures contained entirely within a geological layer.
  • the one or more microfractures are circular and the characterising dimension of each of the microfractures is a radius.
  • the dynamic geomechanical model comprises one or more macrofractures, whose upper and lower edges lie on the upper and lower surfaces of a geological layer.
  • the macrofractures are rectangular or polygonal in shape and the characterising dimension of each of the macrofractures is a length in a direction parallel to the layer by the surfaces of which the macrofracture is limited.
  • the implicit fracture modelling procedure comprises the steps of
  • the data describing a geomechanical model comprises the following:
  • the step of calculating total fracture densities is performed by repeating the following steps for each timestep until the fracture network of the gridblock is complete:
  • the estimation of the optimal timestep duration is done by calculating, for each fracture set, the time taken for the macrofracture density to grow by a specified amount (typically between 0.2% and 1 %) if there is no fracture deactivation and then choosing the smallest of the calculated time values.
  • the calculation of the mean fracture driving stress, the growth factor and the maximum macrofracture propagation distance is based on the estimated optimal timestep duration.
  • the calculation of the rate of macrofracture deactivation is based on the calculated growth factor, the calculated maximum macrofracture propagation distance and the total macrofracture populations as calculated at the end of the previous timestep.
  • the calculation of the total increase in active and static macrofracture populations during the timestep is based on the estimated optimal timestep duration, the calculated mean fracture driving stress, the calculated rate of fracture deactivation, the calculated growth factor and the calculated maximum macrofracture propagation distance for the timestep.
  • the step of calculating the cumulative macrofracture density distribution functions for each fracture set is performed by repeating the following steps for each timestep:
  • relevant dynamic fracture data such as mean fracture driving stress, growth factor, maximum macrofracture propagation distance and macrofracture deactivation rate
  • the step of calculating the cumulative microfracture density distribution functions for each fracture set is performed by repeating the following steps for each timestep:
  • relevant dynamic fracture data such as mean fracture driving stress, growth factor and microfracture deactivation rate
  • implicit fracture data and relevant dynamic fracture data such as mean fracture driving stress, growth factor and maximum macrofracture propagation distance
  • microfractures i. if microfractures are to be included in the explicit discrete fracture network:
  • microfractures are not to be included: 1. adding new macrofractures directly
  • the step of adding new microfractures is performed by:
  • a calculating how frequently microfractures will grow to a specified minimum radius, based on implicit fracture data (such as cumulative density distribution functions) and relevant dynamic fracture data (such as mean fracture driving stress and growth factor) for the relevant fracture set and timestep;
  • implicit fracture data such as cumulative density distribution functions
  • relevant dynamic fracture data such as mean fracture driving stress and growth factor
  • the step of calculating microfracture growth is performed by repeating the following steps for all existing and new microfractures: a. checking if the location of the microfracture lies within a stress shadow of a macrofracture of the same fracture set;
  • microfracture If the microfracture is active: i. increment the radius of the microfracture, based on dynamic fracture data (such as mean fracture driving stress and growth factor) for the relevant fracture set and timestep;
  • dynamic fracture data such as mean fracture driving stress and growth factor
  • the step of adding new macrofractures directly is performed by:
  • the step of calculating macrofracture growth is performed by repeating the following steps for both macrofracture tips of all macrofractures in the fracture set:
  • the step of constructing the explicit discrete fracture network comprises the steps of:
  • each macrofracture must only be checked for intersection and stress shadow interaction against other macrofractures within the same gridblock, not against all other macrofractures in the model.
  • it allows the timestep duration to be optimised independently in each gridblock. This makes the method scalable, allowing it to be applied to small models (e.g. near wellbore models), which must run quickly, as well as to large models (e.g. full field models) that must contain large numbers of fractures.
  • the invention in another aspect of the invention, relates to the use of the method as described above to generate an accurate representation of the fracture network in the subsurface, including key geometric aspects (such as fracture densities, cumulative fracture density distribution functions, fracture porosity and fracture connectivity) for an improved simulation of fluid flow through fractured rock, that can be used to manage oil and gas reservoirs, predict groundwater flow and contaminant dispersion for pollution monitoring, nuclear waste storage, predict flow and leakage in reservoirs used for carbon capture and storage (CCS) and to predict water/fluid flow rates in geothermal systems.
  • key geometric aspects such as fracture densities, cumulative fracture density distribution functions, fracture porosity and fracture connectivity
  • CCS carbon capture and storage
  • it relates to the use of the method as described above to facilitate the generation of accurate, risk-weighted subsurface flow models that can replicate existing oil and gas production data (history matched models), by identifying the key mechanical properties and geological parameters controlling fracture development, and generating multiple geologically plausible fracture models.
  • it relates to the use of the method as described above to predict the impact of natural fractures on the bulk rock mechanical properties for use in assessing mine, tunnel or overland construction stability, for hydrocarbon exploration and production or for monitoring production-related subsidence or risk of seismic activity resulting from oil and gas production, CCS or CO 2 injection.
  • memory operative coupled to the processor, one or more modules of which memory comprises stored processor-executable instructions arranged to instruct the computer system to perform the method as described above.
  • memory operative coupled to the processor, one or more modules of which memory comprises stored processor-executable instructions arranged to instruct the computer system to perform the method as described above.
  • it relates to one or more non-transitory computer-readable storage media comprising computer-executable instructions arranged to instruct a computer system to perform the method as described above.
  • Fig. 1 is a schematic illustration of the single layer fracture model showing
  • Fig. 2 is a schematic illustration showing a single set of Mode 1
  • Fig. 3 is a schematic illustration of macrofracture stress shadow interaction
  • Fig. 4 is a schematic illustration of orthogonal macrofracture intersection
  • Fig. 5 is a graphical representation showing the growth of active half- macrofractures through time according to an embodiment of the invention
  • Fig. 6 is a graphical representation for supporting the calculation of the static volumetric half-macrofracture density according to an embodiment of the invention
  • Fig. 7 is a graphical representation for supporting the calculation of the static mean linear half-macrofracture density according to an embodiment of the invention
  • Fig. 8 shows the overall process steps of an implicit fracture modelling procedure according to an embodiment of the invention
  • Fig. 9 shows the process steps to calculate total fracture densities according to an embodiment of the invention.
  • Fig. 10 shows the process steps to calculate cumulative macrofracture density distribution functions according to an embodiment of the invention
  • Fig. 11 shows the process steps to calculate cumulative microfracture density distribution functions according to an embodiment of the invention
  • Fig. 12 shows the overall process steps of an explicit fracture modelling
  • Fig. 13 shows the process steps to add new microfractures according to an embodiment of the invention
  • Fig. 14 shows the process steps to calculate microfracture growth according to an embodiment of the invention
  • Fig. 15 shows the process steps to add new macrofractures directly according to an embodiment of the invention
  • Figs. 16a-b show the process steps to calculate macrofracture growth according to an embodiment of the invention
  • Fig. 17 shows the process steps to construct an explicit discrete fracture
  • Fig. 18 is a schematic illustration of bevelling of corner points between two non-parallel fracture segments according to an embodiment of the invention. DETAILED DESCRIPTION OF THE INVENTION
  • a single flat, homogeneous and isotropic poroelastic layer of uniform thickness h as illustrated in Fig. 1 is
  • Mode 1 vertical dilatant fractures 101 can be either Mode 1 vertical dilatant fractures 101 or Mode 2 inclined shear fractures 102, depending on the effective vertical stress.
  • Mode 2 fractures 102 dip in one direction or the other at a dip angle w compared to the upper and lower horizontal surfaces of the layer.
  • Mode 2 fractures 102 are assumed to dip in either direction with equal probability, so the net strain on the whole fracture set is horizontal pure shear with no component of simple shear.
  • references to the height and area of Mode 2 fractures 102 are to be taken to mean the height and area projected onto a vertical plane parallel to fracture strike.
  • the rock mass is assumed to be a perfectly elastic material, so the horizontal strain can be accommodated either by elastic strain in the rock mass or by elastic displacement on the fractures.
  • the stress and strain tensors are related by the 3-dimensional form of Hooke’s Law.
  • the fracture network will generally not be, although since the two fracture sets strike perpendicular to the principal horizontal strains, the principal stresses and strains will always be coaxial. Therefore, the orthotropic version of Hooke’s Law can be used to relate the stress and strain tensors, when the fractures have a significant impact on the bulk rock elastic properties of the layer.
  • the layer contains an initial seed population of small circular microfractures pF, following a power-law cumulative density distribution function, as shown by the circles in Fig. 2.
  • These microfractures pF will grow in response to the in situ stress in the layer, induced by the applied strain. Initially, they are contained entirely within the layer, and they propagate at an equal rate in all directions, maintaining a circular shape, until they will span the entire layer. The propagation rates are dependent on in situ stress, fracture size and the mechanical properties of the layer, and are calculated based on published laboratory studies of fracture propagation.
  • the microfractures pF will span the layer when their radius r reaches half the layer thickness h/2, although in practice they may be centred anywhere in the layer, so will often not intersect the top and bottom of the layer simultaneously.
  • Macrofractures MF are assumed to be rectangular with a uniform and constant height equal to the layer thickness h, as shown by the rectangles in Fig. 2. Although they cannot propagate vertically, they can propagate laterally (parallel to the layer) in response to the in situ stress, so they will have a variable length L.
  • Each macrofracture MF will have two lateral tips LT, representing the two opposite ends of the macrofracture MF moving parallel to the layer. The two lateral tips LT of each macrofracture MF propagate independently of each other, and each will cease propagating if they fall into the stress shadow of another parallel macrofracture MF (as indicated by letter A in Fig.
  • the rock mass can also deform by time-dependent inelastic processes such as grain sliding, creep and pressure solution, so that viscoelastic model must be used to relate the stress and strain tensors.
  • the boundary conditions are defined in terms of horizontal applied strain rates ⁇ Gp ⁇ h and ⁇ p Ghqc , rather than instantaneous horizontal strains Shmin and Chmax.
  • the evolution of the fracture network is discretised into a series of timesteps, recalculating the in situ stress, the fracture density and the bulk rock elastic properties at the start of each timestep.
  • Layers can also be modelled with more complex geometries (for example dipping layers of variable thickness) by discretising the layer laterally into a series of gridblocks, as in a static geomodel. If the layer has a relatively shallow dip ( ⁇ c.20°) and its depth and thickness do not vary significantly across a gridblock, reasonable results can be obtained by treating the gridblock as if it were horizontal and planar, using the average depth and thickness to calculate the fracture growth, and then projecting the resulting fracture network onto the true top and bottom surfaces of the gridblock. However, it should be ensured that the length and width of the gridblocks are significantly greater than their height, and (ideally) greater than the typical length of the layer-bound macrofractures MF.
  • More steeply-dipping layers can still be modelled, as long as it is assumed that the fracture network developed when they were relatively flat-lying and they were only subsequently tilted. This is often likely to be a reasonable assumption, as fracture networks often develop early in the deformation episode at relatively low strains.
  • the gridblock must first be reoriented back to horizontal, and the fracture growth must be calculated in a horizontal layer, before the resulting fracture network is rotated back to its current orientation.
  • Microfractures pF are defined as small circular fractures that are contained within the layer. They propagate at an equal rate in all directions, thus remaining circular, until they reach the top and the bottom of the layer. It is assumes that this happens when the fracture radius r reaches half the layer thickness h/2, although in practice the microfractures pF may be centred anywhere in the layer and will often not intersect the top and bottom of the layer simultaneously.
  • a and K c are material properties of the layer: A is the rupture velocity in the material (typically similar to the sonic velocity, e.g. c.2000m/s for consolidated rock), and K c is the critical stress intensity (also known as the fracture toughness), given by where G c is the crack surface energy of the material (typically c.1000 J/m2), E is the Young’s Modulus (c.10GPa) and v the Poisson’s ratio (c.0.25).
  • b represents the subcritical fracture propagation index, which controls the rate of fracture propagation. For low values of b ( ⁇ c.5), slow subcritical fracture propagation will predominate, whereas at high values of b (>.20), rapid critical fracture propagation will predominate.
  • K f represents the stress intensity at the fracture tip, which will be the circumference of the fracture for circular fractures.
  • the stress intensity factor Ki at the tip of a circular Mode 1 dilatant fracture 101 is derived as [3]
  • the effective normal stress o n ’ acting on the fracture is tensile (i.e. negative).
  • the stress intensity factor for a circular Mode 2 shear fracture 102 is more complex, as it is dependent on the magnitude of the shear stress t acting on the fracture and the frictional traction on the fracture surface.
  • the analytical solutions derived for stress and strain around Mode 1 fractures 101 are also valid for Mode 2 fractures 102, if the shear stress minus the frictional traction is substituted for the tensile normal stress in [3].
  • the frictional traction is given by the effective normal stress o n ’ times the fracture friction coefficient p fr , which is a function of the surface roughness. Therefore, the stress intensity factor K 2 for a circular Mode 2 shear fracture 102 can be expressed as
  • microfracture growth equation [11] can be used to derive expressions for the cumulative microfracture density distribution functions MF P3 O and MF P32.
  • the volumetric microfracture density can be described by a cumulative population distribution function MF P3o(r N ,t N ), which gives the number of microfractures of radius G N or greater per m 3 volume at the current time t N .
  • volumetric microfracture density given by [14] will include those microfractures pF that have reached the layer boundaries and become layer-bound macrofractures MF.
  • layer-bound macrofracture population must be subtracted from [14]
  • the total area of all microfractures in a unit volume, MF P 32 (0,t), can be calculated by integrating the number of microfractures of a given radius r times the area of a microfracture with radius r for radii between 0 and h/2 (the maximum radius of a microfracture in the layer).
  • the total area of microfractures per unit volume is equivalent to the number of microfractures that will be encountered along a linear transect perpendicular to the fractures and, thus, can also be called the mean linear microfracture density.
  • the number of microfractures with radius between r and r+dr can be obtained by differentiating the cumulative volumetric microfracture density function MF P o(r,t) with respect to r. Therefore, the total area of all microfractures pF in a unit volume can be expressed as
  • d MF P o/dr at time t N can be obtained by differentiating [14] with respect to r:
  • the cumulative linear microfracture density as a function of microfracture radius, MF P32(r N ,t N ), can be calculated by adjusting the lower limit of the numerical integration in [18].
  • Modelling macrofractures In the following, it is examined how of circular microfractures pF can develop into layer-bound macrofractures MF, and expressions for the rate of macrofracture nucleation and growth are derived. These are used to calculate the cumulative density distribution functions describing the layer-bound macrofracture population, again ignoring fracture interaction.
  • the fracture will span the entire layer and can no longer propagate vertically. It can, however, continue to propagate laterally, i.e. along the layer. It can thus be approximated as a long rectangular layer-bound fracture, which will be labelled as a macrofracture MF.
  • Macrofractures MF will propagate independently at both lateral tips LT until propagation ceases for some reason (for example they intersect an orthogonal macrofracture MF or they enter the stress shadow of another parallel macrofracture MF).
  • the key parameter characterising macrofracture size is thus the macrofracture length L.
  • macrofractures MF can cease propagating (e.g. if they intersect with an orthogonal macrofracture MF or they enter the stress shadow of another parallel macrofracture MF).
  • volumetric macrofracture density MF P 30 The volumetric macrofracture density can be described by a cumulative population distribution function MF P3o(L N ,t N ), which gives the number of macrofractures of length L N or greater per m 3 volume at the current time t N .
  • the number of macrofractures of length LN or greater at time t N is equal to the total number of macrofractures at time t m , which in turn is equal to the total number of microfractures with radius h/2 or greater at time t m , M FP 3 o(h/2,t m ):
  • timestep M is defined such that
  • the mean linear macrofracture density MF P32 can be calculated in the same way as for microfractures pF, by differentiating the volumetric macrofracture density MF P30 with respect to macrofracture length L to obtain the number of macrofractures MF with length between L and L+dL, multiplying the result by fracture area hl_, and integrating across the full range of possible macrofracture lengths L.
  • the fracture height is constant, so the macrofracture area is linearly proportional to the macrofracture length L. This simplifies the calculation of MF P32, as there is an analytical solution to the resulting differential equation:
  • timestep K represents the timestep of nucleation of a macrofracture with current length LN, defined such that
  • [31] can be simplified by cancelling terms across the timesteps, to obtain
  • [33] can be converted into an expression for the total porosity ⁇ 1>MF of a set of Mode 1 dilatant macrofractures 101 if the mean fracture aperture is known. It has been shown that if a cross section is taken of a laterally extensive elastic Mode 1 fracture 101 of height h (i.e. a 2-dimensional plane strain approximation), it will dilate with an elliptical profile, with aperture such that the total volume of a macrofracture of length L, VMF(L), will be given by
  • the total macrofracture porosity at time t N CP MF ⁇ N
  • the total fracture porosity ⁇ 11 ⁇ 2F for a set of vertical Mode 1 dilatant macrofractures 101 also represents the total horizontal extensional strain accommodated by displacement on the fracture set, MFCh (t).
  • the total horizontal extensional strain accommodated by displacement on the fracture set is represented by the total volumetric heave in a unit volume P MR ( ⁇ N ).
  • Fractures are often surrounded by a stress shadow, which comprises the zone in which a change in the applied strain is accommodated by an elastic displacement on the fracture, rather than by a change in the elastic strain (and hence stress) in the rock mass surrounding the fracture.
  • the total width of the stress shadow is given by the change in fracture displacement divided by the change in bulk rock elastic strain. Therefore, [38] can be combined with an expression for the horizontal stress (i.e. the driving stress) as a function of applied horizontal strain, to obtain an expression for the stress shadow width around a macrofracture.
  • the relationship between driving stress and applied horizontal strain in an elastic rock mass can be obtained by combining [6] with the 3-dimensional form of Hooke’s Law, giving for a Mode 1 macrofracture 101 perpendicular to the minimum horizontal applied strain Shmin, and
  • the mean stress shadow width W around individual macrofractures is thus given by for Mode 1 macrofractures 101 with mean aperture a, and
  • [38] can also be combined with [41] and [42] to obtain expressions for the total macrofracture stress shadow volume IJJ MF as a proportion of total rock volume: for Mode 1 dilatant fractures 101 , and
  • the probability that an active fracture will be deactivated and become a static fracture can be calculated. This may take one of two forms: the instantaneous probability of fracture deactivation represents the probability that a specific fracture will cease propagating within a very short time increment dt, while the cumulative fracture activation probability represents the probability that an initially active fracture will still be active at the end of a longer time interval, e.g. a timestep. If fracture deactivation is a random process, so that the instantaneous probability of fracture deactivation q is independent of time, then the cumulative fracture activation probability Q will follow an exponential decay curve of the form
  • fracture deactivation will be a regular or semi-regular process, where the instantaneous probability of fracture deactivation increases with time. In these cases, it is easier to calculate the probability that a fracture will be deactivated over a longer time interval (i.e. the cumulative activation probability), and then divide this by the time interval to obtain a pseudo-instantaneous or mean probability of fracture deactivation.
  • the volumetric density of active fractures can be calculated by combining the cumulative fracture activation probability with the expressions for M FP 30 or MFP30 derived above. If the resulting active fracture volumetric densities, VP30 and 8 MF P3 O , are then multiplied with the instantaneous fracture deactivation probability (or the mean fracture deactivation probability), the rate at which active fractures cease propagating and become static fractures is obtained. Integrating these deactivation rates through time will therefore give expressions for the volumetric density of static fractures, VP30 and S MF P3 O .
  • the previous procedure of differentiating the P30 functions can be repeated to determine the number of fractures with a given area, multiplying this by the fracture area, and integrating the resultant functions across the range of fracture sizes to obtain functions for the volumetric active and static microfracture ratios and mean linear active and static macrofracture densities, VP32, 8 MFP32, VP32, and S MFP32.
  • the fracture population functions include expressions for the fracture deactivation probabilities but are independent of how these are calculated. This is important because, since fractures cease propagating as a result of interaction with other fractures, the fracture deactivation probabilities will change through time as the fracture populations grow.
  • the timestep approach can be used to include this feedback: the fracture deactivation probabilities at the start of each timestep can be calculated, based on the fracture populations at the end of the previous timestep, and then these can be used to calculate the evolution of the fracture population distribution functions during the subsequent timestep.
  • Microfractures will cease propagating if they fall into the stress shadow of a nearby parallel propagating macrofracture.
  • Macrofractures will also cease propagating if their stress shadow interacts with the stress shadow of another parallel macrofracture.
  • Macrofractures will cease propagating if they intersect an orthogonal or oblique macrofracture.
  • microfracture deactivation probability To simplify the calculation of the microfracture deactivation probability, some assumptions about microfracture behaviour are made: • The microfractures are small compared to the macrofractures, and propagate at a slower rate.
  • microfracture falls within the stress shadow zone of a parallel macrofracture. ⁇ The microfracture will have no effect on the macrofracture, which will
  • microfractures have insufficient size or density to develop a significant volume of stress shadow themselves. Interaction with other microfractures is therefore ignored.
  • the proportion of microfractures that will have been deactivated at any later time t N will be given by the ratio of the macrofracture stress shadow volume to the total volume of the layer at time t N . Therefore, the proportion of microfractures still active at time t N , i.e. the cumulative microfracture activation probability QN, is given by:
  • layer-bound macrofractures may cease propagating due to stress shadow interaction.
  • stress shadow interaction There are several key differences:
  • An exclusion zone can be defined around each macrofracture, equal to twice the stress shadow width W, so that if another macrofracture tip enters that exclusion zone, the two stress shadows will overlap and both fracture tips will be deactivated. Unlike the stress shadows, however, the exclusion zones around the macrofractures can partially overlap. If the macrofractures are randomly distributed, but in such a way that that their stress shadows do not overlap, it can be shown that the total exclusion zone volume X MF will be given by
  • a propagating macrofracture tip LTA of a first macrofracture MFA with a stress shadow 201 will interact with the stress shadow 301 around a second parallel macrofracture MFB with exclusion zone 302 if the macrofracture tip LTA lies within the so-called stress shadow interaction box 303 projected ahead a macrofracture tip LTB of
  • This stress shadow interaction box 303 will have width 2W, where W is the mean stress shadow width, and length LSSIB equal to the macrofracture tip propagation rate (i.e. Od b , cf. [20]) times dt, if the macrofracture tip LTB is static, or twice that if both the macrofracture tips LTA, LTB are propagating towards each other.
  • the probability of stress shadow interaction between the two macrofractures MFA, MFB within time interval dt is therefore given by the volume of the stress shadow interaction box 303 as a proportion of the total “clear zone” volume (i.e. the volume not lying within a macrofracture exclusion zone, given by 1 -XMF).
  • the probability that a propagating macrofracture tip LT will experience stress shadow interaction with any other macrofracture MF within time interval dt will be given by the total volume of the stress shadow interaction boxes 303 of all parallel macrofracture tips LT facing in the opposite direction, whether static or dynamic. This will be proportional to the volumetric macrofracture density MFP3O.
  • a function x( ⁇ N ) is defined, representing the total rate of increase of the combined length of all stress shadow interaction boxes facing in one direction at time t N , weighted to take account of obstruction effects:
  • SM MF P3 O and SU MF P3 O represent the volumetric density of macrofractures that are deactivated due to stress shadow interaction and intersection respectively
  • u V and u V are the modifiers to take account of obstruction effects at the static macrofracture tips.
  • the propagating macrofracture tip LT A will always intersect the orthogonal macrofracture MFc before it reaches the stress shadow 301 of the parallel macrofracture MFB. Therefore, it can approximated that u V « 0 (although it is possible to calculate a more precise value if the macrofractures MF are short).
  • layer-bound macrofractures MF will also cease propagating if they intersect a perpendicular-striking macrofracture MF, as long as the frictional traction on the surface of the intersected fracture is sufficiently low to allow slip. It is assumed that this is always the case, so the probability of macrofracture intersection is determined only by the propagation rate of the intersecting macrofracture tip LT and the population of perpendicular-striking macrofractures MF.
  • the probability of macrofracture intersection is dependent not only on the density of the perpendicular-striking macrofractures MF, but also on their spatial distribution. If the perpendicular macrofractures MF are randomly distributed, fracture intersection will be a random process, governed by an exponential probability function. However, if the perpendicular macrofractures MF are evenly spaced (for example if their distribution is controlled by stress shadows), then fracture intersection will be a regular or semi-regular process. In reality, many fracture sets will be somewhere in between these two end-members: There will be a minimum separation between the macrofractures MF, determined by the stress shadow width, but there will also be larger gaps that are randomly distributed.
  • Macrofracture intersection then becomes a semi-regular process that is especially hard to model, as it will include elements of both regular and random behaviour.
  • a general expression for the probability of macrofracture intersection that is valid for any distribution of macrofractures MF can be derived by considering a macrofracture tip LT from a first set of macrofractures MF propagating towards a second set of perpendicular-striking macrofractures MF.
  • the macrofractures MF of this second set are surrounded by stress shadows of width W, which cannot overlap.
  • W may take any value from 0 (for a fully randomly distributed fracture set) to 1/ J MF P 32 (for an evenly spaced set), where J MF P 32 is the mean linear density of the second set of macrofractures.
  • the propagating fracture tip LT will intersect a macrofracture MF of the second set during At if that macrofracture MF lies within the distance DI that the macrofracture tip LT of the first set will propagate in the time period At.
  • a set of“intersection boxes” with width DI can therefore be defined, extending towards the propagating macrofracture tip LT from each of the perpendicular-striking fractures MF.
  • the probability of fracture intersection during the timestep is thus equal to the total volume of all intersection boxes for the second set of macrofractures MF as a proportion of the total rock volume.
  • the total volume of an intersection box attached to a specific macrofracture MF of the second set with length L is given by LhAI. If the width DI of the intersection boxes is less than the stress shadow width W then none of the intersection boxes can overlap, and the total volume of all the intersection boxes is given by J MF P AI.
  • intersection box volume is greater than the stress shadow width W, those parts of the intersection box that lie further from the macrofracture MF than the stress shadow width W may overlap the intersection boxes around other macrofractures MF.
  • a“one-sided proximity zone” formula must be used to calculate the total intersection box volume, similar to that used to calculate total exclusion zone volume in the previous section.
  • the probability F N that a macrofracture which is active at the start of a timestep N will remain active at the end of the timestep can be found by multiplying the cumulative macrofracture activation functions for stress shadow interaction and for intersection:
  • the instantaneous probability of macrofracture deactivation for the timestep F N is much harder to calculate since, as macrofracture deactivation can occur by a combination of random and regular processes, it will vary throughout the timestep in a complex manner.
  • the mean probability F N of macrofracture deactivation can easily be calculated by dividing the probability of macrofracture deactivation during the timestep by the timestep duration:
  • the cumulative probability FN,M that a macrofracture that is active at the end of a timestep M will remain active at the end of a later timestep N:
  • microfracture deactivation probabilities can be combined with the cumulative microfracture population functions derived above to obtain expressions for the active and static microfracture populations.
  • volumetric active and static microfracture density functions VP30 and VP30 can be calculated by modifying the equation [14] for volumetric microfracture density.
  • the active microfracture density VP30 is obtained by multiplying the total microfracture density by the cumulative microfracture activation probability:
  • the rate of microfracture deactivation is equal to the active microfracture population times the instantaneous microfracture deactivation rate.
  • a constant value for Q must be used throughout each timestep, reflecting the macrofracture population at the end of the previous timestep.
  • the integral must be discretised into a sum of terms representing individual timesteps. Also, the mean probability of microfracture deactivation q « rather than the instantaneous probability of microfracture deactivation q « should be used when calculating the rate of microfracture deactivation during a timestep K. This gives
  • [67] and [68] include microfractures that have reached the layer boundaries and become layer-bound macrofractures.
  • Expressions for the active and static macrofracture populations can be obtained in a similar manner: combining the macrofracture deactivation probabilities with the cumulative microfracture population functions derived above.
  • a half-macrofracture ME represents the portion of a macrofracture MF between the nucleation point and one lateral tip LT; each macrofracture MF therefore comprises two conjoined half- macrofractures ME- Therefore, the cumulative macrofracture population functions derived above can easily be converted into cumulative half-macrofracture population functions, expressed in terms of half-macrofracture length I, by replacing the initial microfracture density coefficient B with 2B and replacing the macrofracture propagation rate dL/dt, given by [20], with the half-macrofracture propagation rate dl/dt, given by
  • Fig. 5 illustrates the growth through time of half-macrofractures ME nucleating at different times, specifically at the end of each timestep (the time tx indicating the end of timestep X, timestep N being the current timestep) as a series of lines on a graph of half-macrofracture length against time.
  • the solid lines represent the growth of half-macrofractures ME that nucleate at the timestep boundaries and remain active until the end of timestep N.
  • macrofractures MF may nucleate at any time, and the dashed line indicates the growth trajectory of a half-macrofracture ME that nucleated in the middle of timestep 7. Note that within each timestep, all half-macrofractures ME grow at the same rate, but that this rate can vary between timesteps (in the example shown in Fig. 5, the rate decreases through time, but this does not have to be the case). Thus, at any time t, the time of nucleation of a half-macrofracture ME of length I can be determined by tracing back along the growth curve.
  • [28] cannot simply be multiplied by a single cumulative fracture activation probability to obtain an expression for the cumulative volumetric density of active half-macrofractures ME. as is the case for microfractures pF. Instead, an expression for the number of currently active half-macrofractures ME as a function of half- macrofracture length I must be derived, and then integrated across the range of possible half-macrofracture lengths I.
  • [28] must first be modified to apply to half-macrofractures ME rather than to macrofractures MF, by substituting I for 2L and 2B for B. Then, the maximum density of half-macrofractures ME with length between I and l+dl, if there is no fracture deactivation, is calculated by differentiating the modified version of [8] with respect to I.
  • the nucleation time t m for a half-macrofracture ME of current length I must also be determined by extrapolating back along the half-macrofracture growth trajectory (the solid lines in Fig. 5).
  • the maximum number of half- macrofractures ME can be multiplied by the cumulative fracture activation probability for a fracture nucleating at time t m , and by the effective microfracture activation probability for the time t m , to obtain the density of active half-macrofractures ME with length between I and l+dl.
  • this calculation can be simplified by discretising the cumulative activation probabilities to the timestep boundaries.
  • FN,M-I is defined as the probability that a half-macrofracture ME that nucleated at the start of timestep M is still active at the end of timestep N
  • QM-I is defined as the probability that a given
  • microfracture pF is active at the start of timestep M.
  • These discretised cumulative activation probabilities are then used to calculate the active densities of all half- macrofractures ME that nucleated at any time during timestep M.
  • the integral equation described above can be replaced with a sum across timesteps, so that [28] becomes
  • timestep J is defined such that
  • Fig. 5 To calculate the cumulative volumetric density of all active half-macrofractures ME with length greater than I N at time t N using continuous fracture activation probabilities, the current active density of all half- macrofractures ME with lengths greater than I N , would be integrated along the vertical arrow at the right side of the diagram. If, however, the fracture activation probabilities are discretised by timestep, this integration is converted into a sum of discrete terms representing each timestep from timestep 1 , when the longest currently active half-macrofractures ME nucleated, to timestep J, when half- macrofractures ME with length I N nucleated.
  • a semi-continuous form of [73] for any time t n within timestep N can also be derived:
  • the active half-macrofracture density times the half-macrofracture deactivation rate is integrated with respect to time, along the horizontal arrow in Fig. 6.
  • the active half-macrofracture density is itself an integral with respect to half-macrofracture length I, so overall it is necessary to integrate across the entire hatched area in Fig. 6.
  • timestep J is such that
  • the active half-macrofracture volumetric density equation [75] is itself an integral with respect to half-macrofracture length I, albeit discretised to a sum of terms representing the half-macrofracture nucleation timesteps.
  • the last term in this sum which represents the nucleation timestep of a half-macrofracture of length IN, varies continuously with time. This is illustrated in Fig.
  • M>J for values of t between -i and tK-(aMFAK,M-i-lN)/( aMFadK b )
  • M J for values of t between t K -(aMFA K,M -i-lN)/( a MF a dK b ) and t K -(aMFA K,M -lN)/( a MF a dK b )
  • M ⁇ J for values of t between tK-(OMFAK,M-lN)/( OMFadK b ) and .
  • MF P32(l,t) and s MF P32(l,t) for active and static half- macrofractures ME can be calculated in the same way as MF P32(L,t): by
  • MF P 32 Active mean linear half-macrofracture density 3 MF P 32
  • the differential of MF P3o(l,t) in [30] can be modified to include the cumulative activation functions FN,M-I and O’M-I , and also to represent half-macrofractures ME instead of macrofractures MF. This gives
  • timestep K represents the timestep of nucleation of a half-macrofracture with current length I N , defined such that This time, however, [88] cannot be simplified by cancelling terms across the timesteps, since the cumulative activation probabilities will be different.
  • the mean linear density for static half-macrofractures ME can be calculated in a similar way: by differentiating the volumetric density functions for static half- macrofractures ME with respect to length I to obtain an expression for the number of half-macrofractures ME with length between I and l+dl, multiplying this by half-length I and height h, and then integrating the resulting functions across a range of half- lengths, up to the maximum possible half-length ax . This is illustrated in Fig. 7.
  • the total value of the integral across all lengths between IN and L ax , and from -i to k, can be calculated by summing the integrals for each of the six types of the ⁇ .[ 3 MFP3O>K,M term, over the intervals in which they are valid (note that some types may not be valid over any interval, for any given combination of K and M).
  • S MF P3 O The functions for S MF P3 O derived above are broken down into multiple terms based on the timestep of half-macrofracture nucleation M and the timestep of deactivation K. Each of these terms can be differentiated and reintegrated independently.
  • the differential of the static half-macrofracture volumetric density with respect to length is and the mean linear density of static half-macrofractures s MF P32(l N , ) is given by
  • the integral in [92] can be broken down into six components, each representing the range of half- macrofracture lengths I covered by a type T term ⁇ [ 3 MFP3O>T, as defined above.
  • 3 MFP32>T can be calculated as follows:
  • [104] can be simplified if the value of IN relative to LK,M-I , LK,M, AK-I,M-I and AK-I,M is known. Again, six different types of terms can be defined, in this case reflecting the value of IN relative to the timesteps M and K (i.e. the lowest term component in the summation, represented by the horizontal arrow in Fig. 7). Type 1 terms
  • an alternative method of calculating the half-macrofracture populations in these final timesteps must be found. This can be done by treating the active fracture population as an equilibrium (or“residual”) population, where the rate of fracture nucleation exactly balances the rate of fracture deactivation.
  • the equilibrium rate of half-macrofracture nucleation can be calculated and used to determine the increase in the static half-macrofracture density during the timestep.
  • the size distribution of the equilibrium half-macrofracture population can also be calculated and used to determine the increase in mean linear half-macrofracture density.
  • This technique is particularly useful when simulating the growth of an anisotropic fracture network, with a secondary fracture set that develops only after the primary fracture set has already stopped growing. In this situation, the deactivation rate of the secondary fracture set will be high throughout its growth, so the residual growth technique is the only technique that can accurately predict its population distribution function.
  • the residual active half-macrofracture population is defined as the population where the rate of half-macrofracture deactivation equals the rate of half-macrofracture nucleation, so that the total number of active half-macrofractures remains constant. Therefore where a rMF P3o(t n ) is the residual active half-macrofracture volumetric density, and MFP 3o(t n ) is the half-macrofracture nucleation rate, at any time t n during a timestep N. Note that the instantaneous probability of macrofracture deactivation F N must be used, rather than the mean probability of macrofracture deactivation FN, as the mean macrofracture lifetime is short compared with the timestep duration.
  • the half-macrofracture nucleation rate MF P30 can be calculated by differentiating [28] with respect to time, then multiplying the resulting equation by the effective cumulative microfracture activation probability O’, and also by a factor of 2 to convert from macrofractures MF to half-macrofractures JV . This gives:
  • the rate of change of Q’ can be calculated from the number of propagating macrofracture tips and the rate of tip propagation: Following a similar derivation to that shown above results in
  • [1 19] is difficult to solve, the problem can be simplified by considering two end- member states: one in which the half-macrofracture nucleation rate is controlled predominantly by changes in the microfracture growth rate, and the other in which the half-macrofracture nucleation rate is controlled predominantly by changes in the effective cumulative microfracture activation probability. It can be determined which of these end-member states to use in any given situation by comparing the residual half-macrofracture volumetric densities calculated for the two end member states, and taking the lowest.
  • [113] can simply be used to calculate the half-macrofracture nucleation rate MF P3 O throughout the timestep. Combining this with [112] gives an expression for the total volumetric density of residual active half-macrofractures at any time t n during timestep N:
  • O’ N should be used rather than O’ N -I in [122]
  • O’ N is dependent on MF P 3 o(0,t N ), so first [122] must be used to calculate MFP3o(0,tN), and then MFP3o(0,tN) must be used to update O’N.
  • the length distribution of the residual active half-macrofracture population can be calculated by inserting the half-macrofracture length I, propagation rate and instantaneous deactivation probability into the exponential fracture activation probability function [45].
  • the number of residual active half- macrofractures of length I N or greater is given by
  • microfracture activation probability Q during the timestep, i.e.
  • G N can be substituted for G N -I in [128], since G is not dependent on MFP 3O(0, ⁇ N) ⁇
  • the active half-macrofracture density is multiplied by the probability of macrofracture deactivation, and the result is integrated throughout the timestep.
  • this calculation is easier for the residual fracture population, since it is assumed that all the active half-macrofractures have nucleated during the current timestep, so it is not necessary to iterate through previous timesteps of fracture nucleation.
  • there will be two forms of the equations for residual volumetric density of static macrofractures one form describing residual fracture populations controlled by microfracture growth rate, and one form describing residual fracture populations controlled by the cumulative microfracture activation probability. The form that gives the smallest increase in fracture density will indicate which of these two factors is actually controlling the rate of fracture nucleation at any time.
  • the mean linear density of residual half-macrofractures can be calculated in the same way as above: by differentiating the appropriate MF P30 functions with respect to I to determine the number of half-macrofractures with length between I and l+dl, multiplying this by Ih to obtain their total area, and then integrating the resulting expression from IN to the maximum possible half-macrofracture half-length ax.
  • this is much simpler for residual half-macrofractures, because the cumulative density distribution functions for residual half-macrofractures all follow the form of [123]:
  • the total mean linear density of all residual active half-macrofractures at the end of timestep N is given by and the mean linear density of residual active half-macrofractures of length IN or greater at the end of timestep N is given by where a rMF P 3 o(0,t N ) is given by [122] or [128] as appropriate.
  • This method takes into account the effect that the growing fracture populations will have on the internal elastic stress field within the layer (the in situ stress). This is important because the in situ stress drives fracture propagation, and controls the rate of fracture propagation through the fracture driving stress O d .
  • the growth of the fractures will also affect the in situ stress because fracture displacement accommodates some of the horizontal strain that would otherwise be accommodated elastically in the rock mass, either locally (in the form of a stress shadow) or distributed throughout the layer. This creates a feedback mechanism that is likely to be an important control on the final fracture density and distribution.
  • the boundary conditions for tectonic deformation can typically be defined in terms of a constant effective vertical stress o v ’, and a bimodal horizontal strain applied at specified rates ⁇ Gh ⁇ h and ⁇ p Ghqc , in a specified orientation and for a specified duration.
  • this information combined with the bulk rock elastic moduli, is sufficient to calculate the driving stress on each fracture set, at any time t. If the horizontal strain is applied at a constant rate, the fracture driving stress will increase linearly with time.
  • time-dependent (viscoelastic) deformation processes such as creep, grain sliding and pressure solution, which gradually convert elastic strain into inelastic strain.
  • O dN -i is the fracture driving stress at the end of the previous timestep N-1
  • do d /dt is the rate of change of fracture driving stress during the current timestep N. Since the rate of fracture propagation is proportional to the fracture driving stress to the power of the subcritical fracture propagation index b, the weighted mean driving stress o dN can be substituted for O dN , to calculate the overall growth of the fracture populations throughout the timestep.
  • the bulk rock elastic moduli take account of the elastic strain accommodated on the fractures as well as the elastic strain in the rock mass. Therefore, they may themselves change through time as the fracture network grows. The exact form of the bulk rock elastic moduli will depend on the distribution of stress throughout the layer: whether there are fracture stress shadows or the stress is evenly distributed.
  • e hmm (t) and S hmax (t) represent the instantaneous elastic horizontal strains at time t, which will be a function of the applied strain rates and the strain relaxation rate.
  • the fractures will introduce an anisotropy in the elastic response of the bulk rock to an applied external strain. If it is assumed that the fracture sets are orthogonal to the principal horizontal stress directions, the orthotropic form of Hooke’s Law can be used to relate bulk rock elastic strain to in situ stress:
  • hmin MF P32 represents the mean linear density of all fractures (active and static) striking perpendicular to the Chmin orientation
  • hma) WP32 represents the mean linear density of all fractures striking perpendicular to the S hmax orientation
  • Mhh, M V h, M hv and M w are a set of Fracture Mode factors defined as follows:
  • the method uses a combined implicit and explicit fracture model.
  • the implicit fracture model is used for simulating the processes of fracture nucleation, propagation and interaction, based on geomechanical principles and representing fracture populations statistically by functions or values.
  • the explicit fracture model is used for simulating fracture propagation, interactions and intersections between different fractures in a Discrete Fracture Network (DFN), in which a fracture is represented as an individual object, based on data from the implicit fracture model, such as fracture nucleation, propagation and growth rates.
  • DFN Discrete Fracture Network
  • This complexity and variability can be incorporated into the present fracture models by running the implicit components of the method independently for each gridblock, using the mechanical properties and deformation history defined for that gridblock. This will allow for an easy comparison of the fracture density and geometry predicted in different parts of the geological structure - e.g. the predicted mean fracture spacing and length in a fold hinge and a fold limb, or in the footwall and hanging wall of a fault can be compared.
  • the gridblock geometry there are some constraints on the gridblock geometry: ⁇ The gridblocks should extend vertically across the entire thickness of the fractured layer - i.e. the thickness of the layer should correspond to the height of one gridblock.
  • the lengths and widths of the gridblocks should be greater than the typical maximum length of the macrofractures, so that the majority of macrofractures that nucleate within a gridblock will remain entirely within that gridblock.
  • the thickness and depth of the layer should not vary by too much within each gridblock - i.e. they should be relatively flat-lying (or at least should have been at the time of deformation), and should not have significant taper.
  • Fig. 8 shows the overall process steps of the implicit fracture modelling procedure according to an embodiment of the invention.
  • This procedure is used for generating implicit fracture population data for each gridblock in the static geomodel, in the form of cumulative fracture density functions MF P3o(r), MF P32(r), MFP3O(I) and MFP32(I).
  • the procedure can also be used to generate other useful data described in the previous sections, such as the fracture aperture and porosity, total fracture strain, mean fracture length, and fracture connectivity. This data can be broken down by fracture set if required.
  • Step 801 The implicit fracture modelling procedure is initiated by retrieving data describing the geomechanical model.
  • data includes, for instance, vertical and lateral dimensions of each of the geological layers, mechanical property data for the geological layers and stress and strain history data.
  • Step 802 The model is discretised into a number of gridblocks.
  • Steps 803 and 808 The remaining steps 804-807 are worked through for all gridblocks, one by one.
  • Step 804 The initial stress state is calculated, based on the depth of burial, on the fluid pressure and on the degree of initial compaction. Equations [141]-[151] can be used for this calculation.
  • Step 805 The evolution of the total active and static fracture densities M FP3O(0), M FP 32(0), MFP3O(0) and MFP32(0) for each fracture set are calculated, looping through the timesteps.
  • Step 806 The cumulative macrofracture populations are calculated, looping through the timesteps.
  • Step 807 The cumulative microfracture populations are calculated, looping through the timesteps. Steps 805-807 are described in more detail in the following in relation to Figs. 9-1 1 , respectively.
  • Fig. 9 shows in more detail the process steps to calculate total fracture densities according to an embodiment of the invention.
  • Fig. 9 is a breakdown of step 805.
  • Steps 901 and 910 The individual steps 902-909 are worked through for all timesteps or until the fracture network is complete, one timestep at a time.
  • Step 902 The stress state is recalculated at the start of each timestep. This recalculation is based on the vertical stress, on the fluid pressure, on the total horizontal strain and on the strain relaxation. It also takes account of the impact of fracture growth on the mechanical properties, using fracture populations from the previous timestep. Like in step 804, in which the initial stress state is calculated, Equations [141]-[151] can be used for this calculation.
  • Step 903 The optimal timestep duration is estimated for each timestep.
  • the duration of each timestep is determined dynamically to ensure that the maximum increase in macrofracture stress shadow volume remains below a specified limit; this ensures that the used approximation that fracture deactivation rates are constant throughout the timestep is valid.
  • the estimation is can be done by calculating, for each fracture set, the time taken for the macrofracture density to grow by a specified amount (typically between 0.2% and 1 %) if there is no fracture deactivation and then choosing the smallest of the calculated time values.
  • Step 904 Based on the recalculated stress state, on the estimated optimal timestep duration, on the applied strain rate and on the fracture set orientation, the mean fracture driving stress, i.e. the local stress acting on each fracture set (Equation [140]), the growth factor parameters y and G (Equations [7]-[12]) and the maximum macrofracture propagation distance O MF A (Equations [19]-[29]) are calculated.
  • the mean fracture driving stress i.e. the local stress acting on each fracture set (Equation [140])
  • the growth factor parameters y and G Equations [7]-[12]
  • O MF A Equations [19]-[29]
  • Step 905 The rate of macrofracture deactivation is calculated.
  • Macrofracture deactivation can occur by two mechanisms: Interaction with the stress shadow of other parallel macrofractures or intersection with other non-parallel macrofractures.
  • the probability that a specific macrofracture will be deactivated during a given timestep is a function of the fracture growth rate and the fracture populations (at the end of the previous timestep). Equations [50]-[65] can be used for the calculation of the macrofracture deactivation probabilities F and F.
  • Step 906 The increase in the total active and static macrofracture densities MF P 3O (0) and MF P 32(0) are calculated analytically for each fracture set. These calculations are based on the timestep duration, on the driving stress, on the growth factor, and on the macrofracture deactivation probability for each fracture set. Equations [72]-[11 1] can be used for these calculations. If the macrofracture deactivation probability is high, it may be necessary to use the formulae for residual macrofracture populations (Equations [1 12]-[139]).
  • Step 907 The rate of microfracture deactivation is calculated.
  • Microfracture deactivation can occur when a microfracture falls into the stress shadow of a macrofracture. Thus, it is a function of the growth of the macrofracture population in the timestep and it can only be calculated after the calculation of the macrofracture growth. Equations [45]-[49] can be used for the calculation of the microfracture deactivation probabilities q and Q.
  • Step 908 The increase in the total active and static microfracture densities mr R 3 o(0) and mr R 32(0) are calculated analytically and numerically, respectively, for each fracture set. These calculations are based on the timestep duration, on the driving stress, on the growth factor, and on the microfracture deactivation probability for each fracture set. Equations [66]-[71] can be used for these calculations.
  • Step 909 The calculation continues until the fracture network is complete, i.e. until all fracture sets have stopped growing or the specified maximum time or timestep limit has been reached. A fracture set is assumed to have stopped growing, if the density MF P3 O (0) is below a specified value, or if the total exclusion zone volume approaches 100%.
  • Fig. 10 shows in more detail the process steps to calculate cumulative
  • Fig. 10 is a breakdown of step 806.
  • Steps 1001 and 1004 The individual steps 1002 and 1003 are worked through for all timesteps, one by one.
  • Step 1002 Dynamic fracture data are retrieved for each timestep. These data include the timestep duration, the fracture driving stress, the growth factor and the macrofracture deactivation probability for each fracture set, all of which were calculated and saved in step 805.
  • Step 1003 The increase in the cumulative active and passive macrofracture densities as a function of length MF P3 O (I) and MF P32(I) are calculated analytically for a range of specified half-lengths I. These calculations are based on the timestep duration, on the driving stress, on the growth factor, and on the macrofracture deactivation probability for each fracture set. Equations [72]-[11 1] can be used for these calculations. If the macrofracture deactivation probability is high, it may be necessary to use the formulae for residual macrofracture populations (Equations [1 12]-[139]). The combination of results for all specified half-macrofracture lengths I define the cumulative macrofracture length distribution functions.
  • Fig. 1 1 shows in more detail the process steps to calculate cumulative microfracture densities according to an embodiment of the invention.
  • Fig. 1 1 is a breakdown of step 807.
  • Steps 1101 and 1104 The individual steps 1102 and 1 103 are worked through for all timesteps, one by one.
  • Step 1 102 Dynamic fracture data are retrieved for each timestep. These data include the timestep duration, the fracture driving stress, the growth factor and the microfracture deactivation probability for each fracture set, all of which were calculated and saved in step 805.
  • Step 1 103 The increase in the cumulative active and passive microfracture densities as a function of radius u F P 3 o(r) and ⁇ 32 ( are calculated for a range of specified fracture radii r. These calculations are based on the timestep duration, on the driving stress, on the growth factor, and on the microfracture deactivation probability for each fracture set. ⁇ 30 ( is calculated analytically and ⁇ 32 ( is calculated by numerical integration between the specified radii r. Equations [66]-[71] can be used for these calculations. The combination of results for all specified radii r define the cumulative microfracture size distribution functions.
  • Fig. 12 shows the overall process steps of the explicit fracture modelling procedure according to an embodiment of the invention. This procedure is used for generating the necessary data for construction of the explicit DFN model, based on the calculated implicit fracture data calculated for each of the gridblocks in the geomodel as described above.
  • the DFN model can be populated either with both
  • microfractures and macrofractures or with macrofractures only are added to the DFN according to the macrofracture nucleation rates calculated with the implicit fracture data. This will reduce the total number of fractures in the DFN while retaining all the long fractures that are likely to act as major fluid flow pathways, thus making it more suitable as input for fluid flow modelling.
  • the procedure uses the timesteps defined for each gridblock when calculating the implicit fracture data. It loops through all of the timesteps in all of the gridblocks in the geomodel in chronological order. Thus, it can use the driving stress, nucleation rate and propagation rate data calculated for each fracture set in each timestep in the implicit calculation phase to control the nucleation and growth of the fractures in the explicit DFN.
  • Step 1201 The explicit fracture modelling procedure is initiated by retrieving implicit fracture data and relevant dynamic fracture data for each fracture set and each timestep from the calculated implicit fracture model.
  • the dynamic fracture data includes, for instance, mean fracture growing stress, growth factor and maximum fracture propagation distance.
  • Step 1202 A sorted list of calculation elements in the form of“gridblock-timesteps” is created. Each such gridblock-timestep consists of a specific timestep of a specific gridblock. Since the timestep durations are calculated individually for each gridblock, they will vary between gridblocks. Therefore, a list of all such gridblock-timesteps is created and sorted in chronological order by the end times of each gridblock- timestep.
  • Steps 1203 and 1211 The remaining steps 1204-1210 are worked through for all gridblock-timesteps, one by one.
  • Steps 1204 and 1210 The remaining steps 1205-1209 are worked through for all fracture sets, one by one.
  • Step 1205 If microfractures are to be included in the explicit DFN, steps 1206 and 1207 are performed and step 1208 is skipped. If not, steps 1206 and 1207 are skipped and step 1208 is performed.
  • Step 1206 New microfractures are added.
  • Step 1207 The growth of each microfracture (newly added ones as well as already existing ones) is calculated, including a possible transformation into a macrofracture.
  • Step 1208 New macrofractures are added directly into the model.
  • Step 1209 The growth of each macrofracture is calculated, taking into account possible interactions with parallel macrofractures and intersections with non-parallel macrofractures as well as possible crossings of gridblock boundaries.
  • Steps 1206-1209 are described in more detail in the following in relation to Figs. 13- 16, respectively.
  • Fig. 13 shows in more detail the process steps to add new microfractures according to an embodiment of the invention.
  • Fig. 13 is a breakdown of step 1206.
  • Step 1301 It is calculated how frequently microfractures will grow to a specified minimum radius. Since the total microfracture density is infinite for a power law distribution, a minimum radius limit for the microfractures in the DFN must be specified. Equation [14] can be used to determine how many microfractures will reach the minimum radius during a given timestep. This calculation is based on the microfracture growth rate, on the cumulative microfracture density distribution functions and on the gridblock volume.
  • Step 1302 The current time is set to the start of the timestep being processed. If the stress is not constant, a weighted time can be used.
  • Step 1303 A random location for a new microfracture is picked within the gridblock.
  • Step 1304 It is checked whether the picked location lies within a stress shadow of a macrofracture of the same fracture set.
  • Step 1305 If the picked location lies within such a stress shadow, no new microfracture will be added to the model.
  • Step 1306 If the picked location does not lie within such a stress shadow, a new microfracture is created at the picked location.
  • the nucleation time of the new microfracture is set to be the current time (or alternatively the weighted time) and the radius of the new microfracture is set to be the specified minimum radius.
  • Step 1307 The current time is incremented according to the calculated frequency at which microfractures will grow to a specified minimum radius. Alternatively, the weighted time is used.
  • Step 1308 Steps 1303-1307 are repeated until the current time (or alternatively the weighted time) has reached the end of the timestep being processed.
  • Fig. 14 shows in more detail the process steps to calculate microfracture growth according to an embodiment of the invention.
  • Fig. 14 is a breakdown of step 1207.
  • Steps 1401 and 1410 The individual steps 1402-1409 are worked through for all microfractures, newly added as well as already existing, one by one.
  • Step 1402 It is checked whether the microfracture lies within a stress shadow of a macrofracture of the same fracture set. The microfracture can have fallen into such a stress shadow since the previous timestep.
  • Step 1403 If the microfracture lies within such a stress shadow, the microfracture is deactivated.
  • Step 1404 If the microfracture is not active, steps 1405-1409 are skipped
  • Step 1405 The microfracture radius is incremented.
  • the new microfracture radius is a function of the current radius, of the microfracture growth rate and of the growth time, the growth rate being determined by equation [8] and the fracture driving stress for the timestep being processed.
  • the growth time is the timestep duration, whereas, for microfractures added during the timestep (in step 1306), the growth time is the time passed from the nucleation time to the end of the timestep.
  • Step 1406 It is checked whether the microfracture has grown into a macrofracture, i.e. if the new radius is greater than half the layer thickness.
  • Step 1407 If the microfracture has not grown into a macrofracture, steps 1408 and 1409 are skipped.
  • Step 1408 A new macrofracture of zero length is generated at the location of the microfracture.
  • Step 1409 The microfracture is deactivated.
  • Fig. 15 shows in more detail the process steps to add new macrofractures directly according to an embodiment of the invention. Thus, Fig. 15 is a breakdown of step 1208.
  • Step 1501 It is calculated how frequently macrofractures will nucleate from growing microfractures, based on data from the implicit fracture model. New macrofractures will nucleate whenever a microfracture grows to a radius of half the layer thickness, cf. step 1406 above.
  • Step 1502 The current time is set to the start of the timestep being processed. If the stress is not constant, a weighted time can be used.
  • Step 1503 A random location for a new macrofracture is picked within the gridblock.
  • Step 1504 It is checked whether the picked location lies within an exclusion zone of a macrofracture of the same fracture set.
  • Step 1505 If the picked location lies within such an exclusion zone, no new macrofracture will be added to the model.
  • Step 1506 If the picked location does not lie within such an exclusion zone, a new macrofracture is created at the picked location.
  • the nucleation time of the new macrofracture is set to be the current time (or alternatively the weighted time) and the length of the new macrofracture is set to be zero.
  • Step 1507 The current time is incremented according to the calculated frequency at which macrofractures will nucleate. Alternatively, the weighted time is used.
  • Step 1508 Steps 1503-1507 are repeated until the current time (or alternatively the weighted time) has reached the end of the timestep being processed.
  • Figs. 16a and 16b show in more detail the process steps to calculate macrofracture growth according to an embodiment of the invention. Thus, Figs. 16a and 16b are a breakdown of step 1209.
  • Steps 1601 and 1619 The individual steps 1602-1618 are worked through for all macrofracture tips, of newly added macrofractures as well as of already existing macrofractures, one by one.
  • Step 1602 If the macrofracture tip is not active, steps 1603-1618 are skipped.
  • Step 1603 The maximum propagation distance of the macrofracture tip is calculated. All macrofracture tips within a specific fracture set propagate at the same rate.
  • the propagation rate can be calculated using equation [20], substituting e for L/2.
  • the maximum propagation distance for a given macrofracture tip within a given timestep is the propagation rate times the propagation time within the timestep.
  • the propagation time is the timestep duration, whereas, for macrofractures added during the timestep (in step 1506), the propagation time is the time passed from the nucleation time to the end of the timestep.
  • Step 1604 If the propagating macrofracture tip will interact with a stress shadow of another macrofracture of the same fracture set if growing by the maximum propagation distance, steps 1605-1607 are performed. If not, these three steps are skipped. Such interaction occurs if the stress shadows of the two macrofractures overlap.
  • Step 1605 The maximum propagation distance is reduced to the distance, at which stress shadow interaction occurs.
  • Step 1606 Both of the interacting macrofracture tips are deactivated.
  • Step 1607 If required, a relay segment is added to link the two interacting macrofractures.
  • Step 1608 If the propagating macrofracture tip will intersect a macrofracture of another fracture set if growing by the maximum propagation distance, steps 1609 and 1610 are performed. If not, these two steps are skipped.
  • Step 1609 The maximum propagation distance is reduced to the distance, at which intersection occurs.
  • Step 1610 The propagating macrofracture tip is deactivated.
  • Step 161 1 If the propagating macrofracture tip will cross a gridblock boundary if growing by the maximum propagation distance, steps 1612-1617 are performed and step 1618 is skipped. If not, steps 1612-1617 are skipped and step 1618 is performed.
  • Step 1612 The real time, at which the propagating macrofracture tip will cross the gridblock boundary, is calculated.
  • Step 1613 The maximum propagation distance is reduced to the distance, at which the gridblock boundary is reached.
  • Step 1614 The propagating macrofracture tip is extended by the maximum propagation distance.
  • Step 1615 The propagating macrofracture tip is deactivated.
  • Step 1616 A new macrofracture of zero length is added in the adjacent gridblock at the entering position.
  • the nucleation time is set to the time calculated in step 1612.
  • One of the macrofracture tips of the new macrofracture is connected to the deactivated macrofracture tip at the gridblock boundary.
  • the status of the other macrofracture tip of the new macrofracture is set to active.
  • Step 1617 If fracture propagation has already been calculated up to this time in the adjacent gridblock, the active macrofracture tip of the new macrofracture is extended in the appropriate direction of propagation into the adjacent gridblock until the position corresponding to the latest calculation time of the adjacent gridblock, checking for stress shadow interaction and intersection as described above in steps 1604-1610. It should be noted that the orientation and propagation rate of macrofractures as well as the timestep durations may be different in the adjacent gridblock. Note that if the fracture orientations in the two adjacent gridblocks are convergent, it will be necessary to extend the fracture along the gridblock boundary to avoid crossing backwards and forwards across the gridblock boundary in an infinite loop. If fracture propagation has not yet been calculated up to this time in the adjacent gridblock, the new macrofracture is left at zero length to be calculated in the next timestep for the adjacent gridblock.
  • Step 1618 The propagating macrofracture tip is extended by the maximum propagation distance.
  • Fig. 17 shows the process steps to construct an explicit DFN according to an embodiment of the invention.
  • the macrofractures generated by the above-described modelling procedures are segmented. This allows the DFN to replicate complex macrofracture geometries, for example macrofractures with variable height or curved macrofractures that follow local variations in the in situ stress orientation (as is often observed in outcrop).
  • the connected planar segments must therefore be recombined into single, generally non-planar entities representing entire macrofractures.
  • Step 1701 All macrofractures, which are connected across gridblock boundaries are combined to create continuous macrofractures that span across multiple gridblocks.
  • Step 1702 The positions of the corner points of all macrofractures are calculated to ensure accurate bevelling across gridblock boundaries and macrofracture intersections. Such bevelling is needed for ensuring a smooth alignment of the connected macrofracture segments.
  • Fig. 18 is a schematic illustration of bevelling of corner points between two non- parallel macrofracture segments according to an embodiment of the invention.
  • P1 and P2 are the two top corners of a first macrofracture segment P propagating towards a gridblock boundary GB in a direction from P1 to P2.
  • Q1 and Q2 are the two top corners of a second macrofracture segment Q propagating away from the gridblock boundary in a direction from Q1 to Q2.
  • variables are applied to constrain the scope of the variables (for example, so they apply only to a specific fracture type or set, or a specific timestep). Where these modifiers are absent, the variables can be assumed to apply generally or in situations in which the scope is not defined; e.g. if one fracture type or set is considered, if the timestep is undefined.

Abstract

A method is disclosed for modelling and simulating a fractured geological structure using a dynamic geomechanical model, which is discretised into geological layers, timesteps and gridblocksand comprises an implicit fracture model and an explicit discrete fracture network model, wherein the method comprises an implicit fracture modelling procedure for simulating the processes of fracture nucleation, propagation and interaction, based on fundamental geomechanical principles, to generate an implicit fracture model, in which a fracture population is represented statistically by functions or values, as well as an explicit fracture modelling procedure for simulating the process of fracture propagation in an explicit discrete fracture network model, in which a fracture is represented as an individual object, based on data from the implicit fracture model relating to fracture nucleation and propagation. Furthermore, some uses of such a method and a system performing such a method are disclosed.

Description

A METHOD AND A SYSTEM FOR MODELLING AND SIMULATING A
FRACTURED GEOLOGICAL STRUCTURE
TECHNICAL FIELD
The present invention relates to a method and a system for modelling and simulating a fractured geological structure using a dynamic geomechanical model comprising an implicit fracture model and an explicit discrete fracture network model.
BACKGROUND OF THE INVENTION
The use of models and simulators to improve the understanding, planning and decision making processes with regards to geological structures has been known for many years.
Typically, standard practice in many industries is to develop a computer generated static geological model of a geological structure of interest that can be used to run fluid flow simulations and plan optimal strategies. Such industries include amongst others those working within the fields of oil and gas, geothermal energy,
groundwater, tunnelling, gas sequestration, construction, geological disposal and mining. The geological and physical principles for the modelling and simulation of the geology are the same for all industries, although there may be some differences in the outputs desired depending on the decision to be taken or the action to be performed. At present, such geological models are especially prevalent in the oil and gas industry where they can be used to plan optimal hydrocarbon extraction strategies.
The geological structures of interest often include fractured reservoirs. If there are interbedded layers in the geological structure with contrasting mechanical properties, the fractures may often be confined within brittle layers sandwiched between more ductile layers, forming“layer-bound fractures” and the overall fracture population within each brittle layer can typically be subdivided into sets of approximately parallel fractures. Such fracture networks will have a major impact on fluid flow through a fractured layer, especially when the rock mass has low permeability, as they will provide critical long-distance fluid conduits. They will also have a critical impact on the strength and elastic properties of the fractured layer as a whole (the“bulk rock” properties), potentially introducing planes of weakness and elastic anisotropy to the layer. It is therefore important to incorporate the fracture network into any 3D models of the subsurface geology that may be used for fluid flow modelling, geomechanical modelling or assessing rock stability. Unfortunately, however, most subsurface fractures cannot be imaged or mapped directly as they are below the resolution of geophysical measurements.
The conventional solution to this problem is to try to characterise the subsurface fracture network as best can be done from available data (typically borehole data) and then extrapolate this to build a stochastic fracture model. In the simplest stochastic models (“implicit” fracture models), the fractures are not represented individually, but rather the fracture network is described statistically by one or more spatially variable properties, such as fracture density, fracture orientation, or fracture connectivity. These properties are then used to directly modify key bulk rock properties, such as porosity and permeability (for flow modelling) or failure strength and elastic moduli (for geomechanical modelling). In the case of permeability, a multiplier, known as the Transmissibility Multiplier (TM), for the transmissibility of each gridblock (model cell representing an assumed homogenous physical volume of the modelled geological structure layer) in the static geomodel can be defined, based on the estimated fracture density in that cell. This is the so-called“TM” approach.
The properties that can be measured and extrapolated in this way are usually very limited. For instance, the boreholes provide no data on the fracture sizes or connectivity, so it is not possible to determine the length distribution (or even the mean lengths) of the fractures, nor is it possible to determine how many fractures are interconnected and how many are isolated. Indeed, it is not even possible to determine the average number of fractures per unit volume (the“volumetric fracture density”, known as P30). These properties are all important in controlling the fracture permeability.
A more sophisticated approach, known as the“DFN” approach, is therefore to build an“explicit” stochastic model, comprising a series of geometric objects representing individual fractures, also known as a Discrete Fracture Network (DFN) model. Of course, the fractures in a DFN cannot correspond directly with actual fractures in the subsurface, as these cannot be imaged. Current standard practice is to populate the DFN models stochastically by placing fractures at random locations and with arbitrary sizes until the density of fractures matches an inferred value estimated from well data or otherwise. The fractures are often placed and sized in ways that do not obey geological theory or physical laws.
The aim is to generate a DFN, which has the same overall statistical properties as the real fracture network. Key properties to be replicated are the fracture density, (both the volumetric fracture density (P30) and the mean linear fracture density (P32), the fracture size distribution, and the fracture connectivity.
Unfortunately, neither the fracture size distribution nor the fracture connectivity can be measured directly in the subsurface, even from seismic anisotropy data. It is therefore necessary to rely on measurements from outcrop analogues to condition the stochastic DFN. However, both are emergent properties that derive from the growth of the fracture network. Therefore, a DFN generated using a simple kinematic growth model, i.e. with fractures nucleating and growing at an arbitrary, specified rate and interacting following simple rules, should exhibit a more realistic fracture size distribution and connectivity than a DFN generated by simply placing fractures of arbitrary size in arbitrary locations. It should be remembered, however, that such models are not dynamic, i.e. they do not take into account the stress, applied strain or mechanical properties, and they do not calculate the fracture size distribution or connectivity directly. Rather, these properties must be determined statistically from the explicit DFN.
A third approach, called the mechanical approach, is to generate a fracture model by simulating the growth of the fracture network dynamically, based on fundamental geomechanical principles. The main advantage of this approach is that no a priori knowledge of the fracture density, size distribution, connectivity or any of the other statistical parameters that are so difficult to measure in the subsurface is required. Nor is it necessary to specify the fracture nucleation and propagation rates. All that is needed is a knowledge of the geomechanical properties and deformation history of the fractured layer. This information is usually readily available: the
geomechanical properties can be taken from mechanical tests on core samples, wireline log geophysical data or published data for similar lithologies, while the deformation history is often described in studies of the regional geology, or can be inferred from the geometry of larger-scale geological structures (e.g. major faults, folds or diapirs). DFNs generated in this way will be more realistic than stochastic DFNs as they will by definition be consistent with the geomechanical processes of fracture development, and will hence honour the geology of the fractured layer and the larger-scale geological structure. The statistical properties of the fracture network emerge from the simulation, unlike the stochastic modelling approaches where they are specified in advance. These simulations therefore generate new information, unlike the stochastic DFNs which simply reflecting the (usually unreliable) data used to create them.
In most cases, this approach involves using some numerical geomechanical method to simulate the propagation and interaction of individual fractures and thereby generate an explicit DFN. The finite element method is most commonly used, but other methods such as the boundary element method or the discrete element method have also been used for this purpose. This approach has been used to model artificially induced fractures as well as natural fracture networks. Numerical geomechanical methods allow detailed modelling of local stress anomalies that build up around fractures and faults, as well as local variation in mechanical properties (e.g. cemented nodules, stiffer or more ductile layers, or previous fractures). Thus, they can predict the effect that these local variations will have on the ultimate fracture geometry. However such methods are computationally expensive, as they must include the local stress field around all of the fractures, so it is rarely practical to use them to model more than a few fractures in a restricted area. Certainly it is not possible to generate full DFNs for large geological structures, which may contain hundreds of thousands of individual fractures, using numerical geomechanical methods.
It has therefore been attempted to use simplified (but still geomechanically consistent) rules governing fracture propagation and interaction to generate larger DFNs. Such methods, however, do not directly calculate the cumulative fracture density distribution functions. Instead, these must be derived by a statistical analysis of the explicit DFNs generated by the method. Therefore, in order to build an implicit fracture model using the such methods, it is necessary to first build an explicit DFN, which includes all fractures within the size range of interest. If this step is to be avoided so that the implicit fracture model is generated directly, the mechanical laws governing fracture propagation to the equations describing fracture populations must be applied as a whole (i.e. the cumulative fracture density distribution functions), rather than to individual fractures. This is much more complex as it involves combining statistics and mechanics, but would allow large fracture populations to be modelled very quickly, since the computation time for such calculations will not be proportional to the number of fractures. This technique could therefore be used, for example, to predict the variations in size distribution and connectivity of small, cm- scale fractures across a large geological structure such as a thrust sheet or a salt dome. However, only very few significant developments in this technique are known so far.
Hence, in low permeability geological structures or structures, in which the presence of fractures dominate the flow patterns, the current geological models are
inadequate and inaccurate, and it is well known in the industries that the margin of error is high and that the value of the predictive results is limited.
The same issues are prevalent within a number of industries that require a thorough understanding of fluid flow through one or more geological layers. The low accuracy of current models used in geothermal energy results in an inaccurate determination of the potential energy output from a proposed well. This could result in a
geothermal well being drilled that is not profitable or a geothermal project not proceeding due to an underestimation of the potential of the project. In the gas sequestration industry, for example Carbon Capture and Storage (CCS), it is also important to be able to predict the success of a project to ensure that it is profitable or that opportunities are not missed. In the case of carbon dioxide injection, it is very important to ensure that any potential leak of the C02 gas will be constrained within the geological structure. In some types of projects relating, for instance, to CCS or geological disposal of nuclear waste, a high degree of certainty is crucial, and any potential existence of fractures and potential leak paths will end the project immediately.
BRIEF DESCRIPTION OF THE INVENTION
It is an object of the present invention to provide a method and a system for modelling and simulating a fractured geological structure that overcome at the above-mentioned disadvantages of the approaches known in the art.
The present invention relates to a method for modelling and simulating a fractured geological structure using a dynamic geomechanical model, which is discretised into geological layers, timesteps and gridblocks and comprises an implicit fracture model and an explicit discrete fracture network model, wherein the method comprises an implicit fracture modelling procedure for simulating the processes of fracture nucleation, propagation and interaction, based on fundamental geomechanical principles, to generate an implicit fracture model, in which model a fracture population is represented statistically by functions or values, and an explicit fracture modelling procedure for simulating the process of fracture propagation, interactions between parallel fractures and intersections of non-parallel fractures in an explicit discrete fracture network model, in which a fracture is represented as an individual object, based on data from the implicit fracture model relating to fracture nucleation and propagation (such as rates of fracture nucleation and growth), as well as incorporating the explicit discrete fracture network model into a geological simulator.
In comparison to methods previously known within this field of art, the method according to the present invention does not only give a more realistic representation of the fractures, honouring the geomechanics and the geology but also improves the accuracy of the resulting dynamic geomechanical model of the geological structure of interest and reduces the time taken to generate it.
Thus, the use of a combined implicit and explicit model makes it possible to simulate fracture networks through a full range of different fracture sizes scales and across large geological structures within a reasonable time and using commonly available input data. The method according to the present invention is able to model both Mode 1 dilatant fractures and Mode 2 shear fractures, and to realistically recreate the geometry of the fracture network developed by the interactions between two sets of fractures, striking perpendicular to the minimum and maximum horizontal strain orientations. It is also capable of generating a purely implicit fracture model, in which the fracture sets are described by cumulative fracture density distribution functions, and can be used to generate large fracture models (e.g. across an entire geological structure) very quickly.
The method allows for the determination of the impact of the growing fracture network on the bulk rock elastic properties, and thus the changes in the in situ stress and the rate of fracture nucleation and propagation as the fracture network grows. This dynamic data from the implicit model, combined with the rules governing fracture interaction, can be used to generate a geomechanically-based explicit DFN without simulation of the propagation of each fracture numerically.
Implicit fracture models can be generated independently for every gridblock in a static geomodel, using mechanical property and horizontal strain data specific to that gridblock, and these results can be combined to generate a single, integrated explicit DFN spanning the entire geomodel. This is a key development as it allows for the generation of large, geomechanically consistent DFNs that reflect the variability in geometry, mechanical properties and in situ stress across a large geological structure, for example a DFN containing all layer-bound fractures across an entire salt diapir, or a near wellbore DFN containing the larger microfractures that develop in the damage zone around a nearby fault. This cannot be done using previous techniques known within the art. Other key advantages of this method over previous methods are:
• The combined implicit and explicit geomechanically-based approach allows for the inclusion of fractures through a full range of scales, from mm- or cm- size microfractures to small faults with lengths of hundreds of metres, in a single, geomechanically consistent model.
• Only this combined implicit and explicit geomechanically-based approach can generate geomechanically-based fracture models across large geological structures such as salt domes or large anticlines, as numerical geomechanical methods (even those using simplified mechanical rules) cannot run fast enough to generate such models within a reasonable time.
• The implicit component of this method can calculate cumulative density distribution functions for the volumetric fracture density P3o(S), the mean linear fracture density P32(S), the fracture porosity (P(S) and the total stress shadow volume IJJ(S), for each fracture set, over a full range of scales from mm-scale microfractures to large layer-bound fractures and even small faults. Note that these cumulative density distribution functions will not necessarily follow a power-law relationship, even if the original“seed” microfracture populations do, due to the effects of differential growth rates, fracture interaction and changes in the in situ stress as the fracture network grows.
• Since fracture interactions are included in the implicit component of the method, the fracture connectivity for the implicit fracture populations can be calculated directly, without needing to generate an explicit DFN.
• By generating implicit fracture models independently for every gridblock in the static geomodel, based on mechanical property and in situ stress data specific to that gridblock, and combining the results, a single, integrated explicit DFN can be generated where fracture density and orientation varies laterally, reflecting local variations in the lithology or applied strain. It is even possible to generate a DFN in which the longest layer-bound fractures are curved, following curves in the situ stress orientation around larger faults or other geological structures.
• It is also possible to use this method to study the development of fracture networks through geological time. The growth of individual fractures can be followed, from small circular microfractures to large layer-bound fractures and small faults, and the evolution of the cumulative fracture density distribution functions and fracture connectivity values can be followed through time. Therefore, multiple, geologically-consistent fracture models can be generated quickly, each representing a different stage in the growth of the fracture network through time. These models can be used for uncertainty and risk analysis.
• Finally, the method can be calibrated by comparing the implicit cumulative fracture density distribution functions against the cumulative fracture density distribution functions calculated for the explicit DFN. A close match will demonstrate that the implicit component of the method is correctly predicting the fracture populations, based on the premises concerning fracture nucleation and interaction.
In an embodiment of the invention, the dynamic geomechanical model comprises one or more microfractures contained entirely within a geological layer.
In an embodiment of the invention, the one or more microfractures are circular and the characterising dimension of each of the microfractures is a radius.
In an embodiment of the invention, the dynamic geomechanical model comprises one or more macrofractures, whose upper and lower edges lie on the upper and lower surfaces of a geological layer.
In an embodiment of the invention, the macrofractures are rectangular or polygonal in shape and the characterising dimension of each of the macrofractures is a length in a direction parallel to the layer by the surfaces of which the macrofracture is limited. In an embodiment of the invention, the implicit fracture modelling procedure comprises the steps of
a. retrieving data describing a geomechanical model associated with a geologic environment;
b. discretising the geomechanical model into a number of gridblocks;
c. performing the following steps for each of the gridblocks:
i. calculating an initial stress state;
repeating each of the following until a specified time or a specified number of timesteps is reached:
ii. calculating total fracture densities;
iii. calculating cumulative macrofracture density distribution functions for each fracture set;
iv. calculating cumulative microfracture density distribution functions for each fracture set.
In an embodiment of the invention, the data describing a geomechanical model comprises the following:
a. vertical dimension and lateral dimensions of each of the geological layers; b. mechanical property data for the geological layers;
c. stress and strain history data.
In an embodiment of the invention, the step of calculating total fracture densities is performed by repeating the following steps for each timestep until the fracture network of the gridblock is complete:
a. recalculating the stress state;
b. estimating an optimal timestep duration;
c. calculating and saving a mean fracture driving stress, a growth factor and a maximum macrofracture propagation distance for each fracture set;
d. calculating a rate of macrofracture deactivation for each fracture set;
e. calculating the total increase in active and static macrofracture density
distribution functions for each fracture set;
f. calculating a rate of microfracture deactivation for each fracture set; g. calculating the total increase in active and static microfracture density distribution functions for each fracture set;
h. checking if the fracture network of the gridblock is complete, i.e. if all
macrofractures in all fracture sets have stopped growing or if the specified time or the specified number of timesteps has been reached.
In an embodiment of the invention, the estimation of the optimal timestep duration is done by calculating, for each fracture set, the time taken for the macrofracture density to grow by a specified amount (typically between 0.2% and 1 %) if there is no fracture deactivation and then choosing the smallest of the calculated time values.
In an embodiment of the invention, the calculation of the mean fracture driving stress, the growth factor and the maximum macrofracture propagation distance is based on the estimated optimal timestep duration.
In an embodiment of the invention, the calculation of the rate of macrofracture deactivation is based on the calculated growth factor, the calculated maximum macrofracture propagation distance and the total macrofracture populations as calculated at the end of the previous timestep.
Determining the proportion of macrofractures, which terminate due to the“stress shadow” effect of other parallel macrofractures or due to intersections with other non-parallel macrofractures at any time, results in more accurate estimates of the total density of macrofractures, the spatial distribution of macrofractures, the distribution of different macrofracture sizes and the macrofracture connectivity than in current methods, which ignore this effect.
In an embodiment of the invention, the calculation of the total increase in active and static macrofracture populations during the timestep is based on the estimated optimal timestep duration, the calculated mean fracture driving stress, the calculated rate of fracture deactivation, the calculated growth factor and the calculated maximum macrofracture propagation distance for the timestep. In an embodiment of the invention, the step of calculating the cumulative macrofracture density distribution functions for each fracture set is performed by repeating the following steps for each timestep:
a. retrieving relevant dynamic fracture data (such as mean fracture driving stress, growth factor, maximum macrofracture propagation distance and macrofracture deactivation rate) for the relevant timestep;
b. calculating the increase in the cumulative macrofracture density distribution functions for a range of specified half-macrofracture lengths. In an embodiment of the invention, the step of calculating the cumulative microfracture density distribution functions for each fracture set is performed by repeating the following steps for each timestep:
a. retrieving relevant dynamic fracture data (such as mean fracture driving stress, growth factor and microfracture deactivation rate) for the relevant timestep;
b. calculating the increase in the cumulative microfracture density distribution functions for a range of specified fracture radii.
The method according to any of the preceding claims, wherein the explicit fracture modelling procedure comprises the steps of
a. retrieving implicit fracture data and relevant dynamic fracture data (such as mean fracture driving stress, growth factor and maximum macrofracture propagation distance) for each fracture set and each timestep from the calculated implicit fracture model;
b. creating a sorted list of gridblock-timesteps, each consisting of a specific timestep of a specific gridblock, by generating a list of all timesteps in all gridblock and sorting it by the end times of the timesteps;
c. performing the following steps for each fracture set in each gridblock- timestep:
i. if microfractures are to be included in the explicit discrete fracture network:
1. adding new microfractures;
2. calculating microfracture growth;
ii. if microfractures are not to be included: 1. adding new macrofractures directly
iii. calculating macrofracture growth;
d. constructing the explicit discrete fracture network for the whole grid.
In an embodiment of the invention, the step of adding new microfractures is performed by:
a. calculating how frequently microfractures will grow to a specified minimum radius, based on implicit fracture data (such as cumulative density distribution functions) and relevant dynamic fracture data (such as mean fracture driving stress and growth factor) for the relevant fracture set and timestep;
b. setting current time to start of the timestep;
c. repeating the following steps until the end of the timestep has been reached: i. picking a random location for a new microfracture;
ii. checking if the location for the new microfracture lies within a stress shadow of a macrofracture of the same fracture set;
iii. if the location does not lie within such a stress shadow:
1. adding a new microfracture to the explicit discrete fracture network at the picked location, setting the status of the new microfracture to active;
iv. incrementing the current time by an amount of time corresponding to the calculated frequency of microfractures growing to the specified minimum radius.
In an embodiment of the invention, the step of calculating microfracture growth is performed by repeating the following steps for all existing and new microfractures: a. checking if the location of the microfracture lies within a stress shadow of a macrofracture of the same fracture set;
b. If the location lies within such a stress shadow:
i. deactivating the microfracture;
c. checking if the microfracture is active;
d. If the microfracture is active: i. increment the radius of the microfracture, based on dynamic fracture data (such as mean fracture driving stress and growth factor) for the relevant fracture set and timestep;
ii. checking if the microfracture has grown into a macrofracture;
iii. if the microfracture has grown into a macrofracture:
1. generating a new macrofracture of zero length at the location of the microfracture;
2. deactivating the microfracture.
In an embodiment of the invention, the step of adding new macrofractures directly is performed by:
a. calculating how frequently macrofractures will nucleate from growing
microfractures, based on data from the implicit fracture model;
b. setting current time to start of the timestep;
c. repeating the following steps until the end of the timestep has been reached: i. picking a random location for a new macrofracture;
ii. checking if the location for the new macrofracture lies within a
macrofracture exclusion zone;
iii. if the location does not lie within such an exclusion zone:
1. adding a new macrofracture of zero length to the explicit discrete fracture network at the picked location, setting the status of the macrofracture tips to active;
d. incrementing the current time by an amount of time corresponding to the calculated frequency of macrofractures nucleating from growing
microfractures.
In an embodiment of the invention, the step of calculating macrofracture growth is performed by repeating the following steps for both macrofracture tips of all macrofractures in the fracture set:
a. checking if the macrofracture tip is active;
b. if the macrofracture tip is active:
i. calculating the maximum propagation distance of the macrofracture tip from the dynamic data calculated for the implicit fracture model; ii. checking if the fracture tip will interact with a stress shadow of another macrofracture of the same fracture set if growing by the maximum propagation distance;
iii. if such interaction will occur:
1. reducing the maximum propagation distance to the distance, at which stress shadow interaction occurs;
2. deactivating both macrofracture tips;
3. if required, adding a relay segment to link the two macrofractures.
iv. checking if the fracture tip will intersect a macrofracture of another fracture set if growing by the maximum propagation distance;
v. if such intersection will occur:
1. reducing the maximum propagation distance to the distance, at which intersection occurs;
2. deactivating the propagating macrofracture tip;
vi. checking if the fracture tip will cross a gridblock boundary to an adjacent gridblock if growing by the maximum propagation distance; vii. if such crossing will occur:
1. calculating the real time, at which the macrofracture tip will cross the gridblock boundary;
2. reducing the maximum propagation distance to the distance, at which the gridblock boundary is reached;
3. extending the macrofracture tip by the maximum propagation distance;
4. deactivating the propagating macrofracture tip;
5. adding a new macrofracture of zero length at the entering position within the adjacent gridblock, connecting one of the macrofracture tips to the deactivated macrofracture tip at the gridblock boundary and setting the status of the other macrofracture tip to active;
6. checking if fracture propagation in the adjacent gridblock has already been calculated up to this time;
7. if this is the case: a. extending the active macrofracture tip of the new
macrofracture into the adjacent gridblock in the appropriate direction of propagation until the position corresponding to the latest calculation time of the adjacent gridblock, checking for stress shadow interaction and intersection as described above in steps ii-v;
viii. if no such crossing will occur:
1. extending the macrofracture tip by the maximum propagation distance.
Determining the proportion of macrofractures, which terminate due to the "stress shadow" effect of other parallel macrofractures or due to intersections with other non-parallel macrofractures at any time, results in more accurate estimates of the total density of macrofractures, the spatial distribution of macrofractures, the distribution of different macrofracture sizes and the macrofracture connectivity than in current methods, which ignore this effect.
In an embodiment of the invention, the step of constructing the explicit discrete fracture network comprises the steps of:
a. combining all macrofractures, which are connected across gridblock
boundaries to create continuous macrofractures that span across multiple gridblocks;
b. calculating positions of corner points of all macrofractures to ensure accurate bevelling across gridblock boundaries and macrofracture intersections.
Determining for each gridblock, which individual macrofractures propagate across gridblock boundaries and linking corresponding macrofractures on the sides of gridblock boundaries together to form continuous macrofractures that span across multiple gridblocks significantly reduces the computation time for two reasons:
Firstly, each macrofracture must only be checked for intersection and stress shadow interaction against other macrofractures within the same gridblock, not against all other macrofractures in the model. Secondly, it allows the timestep duration to be optimised independently in each gridblock. This makes the method scalable, allowing it to be applied to small models (e.g. near wellbore models), which must run quickly, as well as to large models (e.g. full field models) that must contain large numbers of fractures.
In another aspect of the invention, it relates to the use of the method as described above to generate an accurate representation of the fracture network in the subsurface, including key geometric aspects (such as fracture densities, cumulative fracture density distribution functions, fracture porosity and fracture connectivity) for an improved simulation of fluid flow through fractured rock, that can be used to manage oil and gas reservoirs, predict groundwater flow and contaminant dispersion for pollution monitoring, nuclear waste storage, predict flow and leakage in reservoirs used for carbon capture and storage (CCS) and to predict water/fluid flow rates in geothermal systems.
In another aspect of the invention, it relates to the use of the method as described above to facilitate the generation of accurate, risk-weighted subsurface flow models that can replicate existing oil and gas production data (history matched models), by identifying the key mechanical properties and geological parameters controlling fracture development, and generating multiple geologically plausible fracture models.
In another aspect of the invention, it relates to the use of the method as described above to predict the impact of natural fractures on the bulk rock mechanical properties for use in assessing mine, tunnel or overland construction stability, for hydrocarbon exploration and production or for monitoring production-related subsidence or risk of seismic activity resulting from oil and gas production, CCS or CO2 injection.
In a further aspect of the invention, it relates to a computer system comprising
a. a processor;
b. memory operative coupled to the processor, one or more modules of which memory comprises stored processor-executable instructions arranged to instruct the computer system to perform the method as described above. In a further aspect of the invention, it relates to one or more non-transitory computer-readable storage media comprising computer-executable instructions arranged to instruct a computer system to perform the method as described above. THE DRAWINGS
In the following a few exemplary embodiments of the invention are described in further detail with reference to the drawings, of which
Fig. 1 is a schematic illustration of the single layer fracture model showing
Mode 1 (vertical) and Mode 2 (inclined) fracture orientations,
Fig. 2 is a schematic illustration showing a single set of Mode 1
microfractures and macrofractures,
Fig. 3 is a schematic illustration of macrofracture stress shadow interaction,
Fig. 4 is a schematic illustration of orthogonal macrofracture intersection, Fig. 5 is a graphical representation showing the growth of active half- macrofractures through time according to an embodiment of the invention,
Fig. 6 is a graphical representation for supporting the calculation of the static volumetric half-macrofracture density according to an embodiment of the invention,
Fig. 7 is a graphical representation for supporting the calculation of the static mean linear half-macrofracture density according to an embodiment of the invention, Fig. 8 shows the overall process steps of an implicit fracture modelling procedure according to an embodiment of the invention,
Fig. 9 shows the process steps to calculate total fracture densities according to an embodiment of the invention,
Fig. 10 shows the process steps to calculate cumulative macrofracture density distribution functions according to an embodiment of the invention, Fig. 11 shows the process steps to calculate cumulative microfracture density distribution functions according to an embodiment of the invention,
Fig. 12 shows the overall process steps of an explicit fracture modelling
procedure according to an embodiment of the invention,
Fig. 13 shows the process steps to add new microfractures according to an embodiment of the invention,
Fig. 14 shows the process steps to calculate microfracture growth according to an embodiment of the invention,
Fig. 15 shows the process steps to add new macrofractures directly according to an embodiment of the invention, Figs. 16a-b show the process steps to calculate macrofracture growth according to an embodiment of the invention,
Fig. 17 shows the process steps to construct an explicit discrete fracture
network according to an embodiment of the invention, and
Fig. 18 is a schematic illustration of bevelling of corner points between two non-parallel fracture segments according to an embodiment of the invention. DETAILED DESCRIPTION OF THE INVENTION
Conceptual model
In the following, the conceptual model of fracture growth, on which the method is based, is described and the simplifications and assumptions it includes are detailed.
In the simplest version of the conceptual model, a single flat, homogeneous and isotropic poroelastic layer of uniform thickness h as illustrated in Fig. 1 is
considered. Fractures are restricted to the layer and cannot propagate out of it. The layer is subjected to a constant and uniform vertical effective stress ov’ and to extensional horizontal principal strains Ohmin and Ohmax (by convention, compressive stresses and strains are positive so the extensional horizontal strains will be negative. Ohmin represents the most extensional, and hence most negative horizontal strain).
Two sets of fractures can develop in response to the applied horizontal strains, striking perpendicular to Ohmin and Ohmax, respectively (see Fig. 1 ). These can be either Mode 1 vertical dilatant fractures 101 or Mode 2 inclined shear fractures 102, depending on the effective vertical stress. Mode 2 fractures 102 dip in one direction or the other at a dip angle w compared to the upper and lower horizontal surfaces of the layer. However, Mode 2 fractures 102 are assumed to dip in either direction with equal probability, so the net strain on the whole fracture set is horizontal pure shear with no component of simple shear. Note also that references to the height and area of Mode 2 fractures 102 are to be taken to mean the height and area projected onto a vertical plane parallel to fracture strike.
The rock mass is assumed to be a perfectly elastic material, so the horizontal strain can be accommodated either by elastic strain in the rock mass or by elastic displacement on the fractures. The stress and strain tensors are related by the 3-dimensional form of Hooke’s Law. However, while the elastic properties of the rock mass are assumed to be isotropic, the fracture network will generally not be, although since the two fracture sets strike perpendicular to the principal horizontal strains, the principal stresses and strains will always be coaxial. Therefore, the orthotropic version of Hooke’s Law can be used to relate the stress and strain tensors, when the fractures have a significant impact on the bulk rock elastic properties of the layer.
It is assumed that the layer contains an initial seed population of small circular microfractures pF, following a power-law cumulative density distribution function, as shown by the circles in Fig. 2. These microfractures pF will grow in response to the in situ stress in the layer, induced by the applied strain. Initially, they are contained entirely within the layer, and they propagate at an equal rate in all directions, maintaining a circular shape, until they will span the entire layer. The propagation rates are dependent on in situ stress, fracture size and the mechanical properties of the layer, and are calculated based on published laboratory studies of fracture propagation. For simplicity it is assumed that the microfractures pF will span the layer when their radius r reaches half the layer thickness h/2, although in practice they may be centred anywhere in the layer, so will often not intersect the top and bottom of the layer simultaneously.
At this point they become layer-bound macrofractures MF. Macrofractures MF are assumed to be rectangular with a uniform and constant height equal to the layer thickness h, as shown by the rectangles in Fig. 2. Although they cannot propagate vertically, they can propagate laterally (parallel to the layer) in response to the in situ stress, so they will have a variable length L. Each macrofracture MF will have two lateral tips LT, representing the two opposite ends of the macrofracture MF moving parallel to the layer. The two lateral tips LT of each macrofracture MF propagate independently of each other, and each will cease propagating if they fall into the stress shadow of another parallel macrofracture MF (as indicated by letter A in Fig. 2), or if they intersect with an orthogonal macrofracture MF (as indicated by letter B in Fig. 2). Eventually, all the macrofractures MF in the layer will cease propagating. At this point, the fracture network will be“fully saturated” and will not develop further, although the macrofractures MF themselves can continue to accommodate strain by accumulating displacement, and the simulation can be ended.
In more sophisticated versions of the model, the rock mass can also deform by time- dependent inelastic processes such as grain sliding, creep and pressure solution, so that viscoelastic model must be used to relate the stress and strain tensors. In this case, the boundary conditions are defined in terms of horizontal applied strain rates έ Gpίh and έpGhqc, rather than instantaneous horizontal strains Shmin and Chmax.
Furthermore, as the fracture network grows, the amount of strain accommodated on the fractures pF, MF will increase. Therefore, the evolution of the fracture network is discretised into a series of timesteps, recalculating the in situ stress, the fracture density and the bulk rock elastic properties at the start of each timestep.
Layers can also be modelled with more complex geometries (for example dipping layers of variable thickness) by discretising the layer laterally into a series of gridblocks, as in a static geomodel. If the layer has a relatively shallow dip (<c.20°) and its depth and thickness do not vary significantly across a gridblock, reasonable results can be obtained by treating the gridblock as if it were horizontal and planar, using the average depth and thickness to calculate the fracture growth, and then projecting the resulting fracture network onto the true top and bottom surfaces of the gridblock. However, it should be ensured that the length and width of the gridblocks are significantly greater than their height, and (ideally) greater than the typical length of the layer-bound macrofractures MF.
More steeply-dipping layers can still be modelled, as long as it is assumed that the fracture network developed when they were relatively flat-lying and they were only subsequently tilted. This is often likely to be a reasonable assumption, as fracture networks often develop early in the deformation episode at relatively low strains. In this case, the gridblock must first be reoriented back to horizontal, and the fracture growth must be calculated in a horizontal layer, before the resulting fracture network is rotated back to its current orientation.
Modelling microfractures
In the following, the growth of circular microfractures pF is considered, and simple expressions for the evolution of the cumulative density distribution functions describing the microfracture population are derived, ignoring fracture interaction. Microfractures pF are defined as small circular fractures that are contained within the layer. They propagate at an equal rate in all directions, thus remaining circular, until they reach the top and the bottom of the layer. It is assumes that this happens when the fracture radius r reaches half the layer thickness h/2, although in practice the microfractures pF may be centred anywhere in the layer and will often not intersect the top and bottom of the layer simultaneously.
Microfracture propagation rate
Experimental studies have shown that the fracture propagation rate dr/dt follows a power law relationship with the stress acting on the microfracture pF, given by
Figure imgf000025_0001
Here, A and Kc are material properties of the layer: A is the rupture velocity in the material (typically similar to the sonic velocity, e.g. c.2000m/s for consolidated rock), and Kc is the critical stress intensity (also known as the fracture toughness), given by
Figure imgf000025_0002
where Gc is the crack surface energy of the material (typically c.1000 J/m2), E is the Young’s Modulus (c.10GPa) and v the Poisson’s ratio (c.0.25). b represents the subcritical fracture propagation index, which controls the rate of fracture propagation. For low values of b (<c.5), slow subcritical fracture propagation will predominate, whereas at high values of b (>.20), rapid critical fracture propagation will predominate.
Kf represents the stress intensity at the fracture tip, which will be the circumference of the fracture for circular fractures. The stress intensity factor Ki at the tip of a circular Mode 1 dilatant fracture 101 is derived as [3]
Figure imgf000026_0001
where the effective normal stress on’ acting on the fracture is tensile (i.e. negative). The stress intensity factor for a circular Mode 2 shear fracture 102 is more complex, as it is dependent on the magnitude of the shear stress t acting on the fracture and the frictional traction on the fracture surface. However, it has been shown that the analytical solutions derived for stress and strain around Mode 1 fractures 101 are also valid for Mode 2 fractures 102, if the shear stress minus the frictional traction is substituted for the tensile normal stress in [3]. The frictional traction is given by the effective normal stress on’ times the fracture friction coefficient pfr, which is a function of the surface roughness. Therefore, the stress intensity factor K2 for a circular Mode 2 shear fracture 102 can be expressed as
Figure imgf000026_0002
Following solutions known within the art, subsequent calculations can be simplified by deriving a generalised expression for the stress intensity factor Kf for any circular microfracture pF:
[5]
Figure imgf000026_0003
where Od is the fracture driving stress, given by for Mode 1 fractures
for Mode 2 fractures
Figure imgf000026_0004
The values of the normal stress on’ acting on the fracture and shear stress magnitude |t| acting on a specific fracture are of course dependent on the fracture orientation and the orientation and magnitude of the in situ stress field. Microfracture growth
Since Kr is a function of the fracture radius r, combining [1] and [5] gives a differential equation
[7] dr ( 2
- dt = A te)
Figure imgf000027_0001
A solution to [7] (when b>2) is given by
[8] r(t) = r(| tj - t\ where
Figure imgf000027_0002
Y is hereafter referred to as the microfracture growth factor, and controls the microfracture growth rate. Note that b is negative when b>2. The form of the solution is slightly different when b=2 or b<2.
[8] can be rearranged to calculate the microfracture radius ro at time t=0 if the radius n at time t=ti is known:
[10] r0 = (ry + gbί1
[10] assumes that the driving stress Od is constant. However, in general Od will vary with time. If time is broken down into a series of discrete timesteps, and it is assumed that Od is constant within each timestep but can vary between timesteps, a form of [10] valid for variable driving stress can be derived: [1 1 ] r0 = (rN p + GN) where GN is the cumulative microfracture growth factor. This will later be used to calculate macrofracture nucleation rate and density. It is defined as
[12]
Figure imgf000028_0001
t -i)] for a series of timesteps J up to the present timestep N, where yj and Odj represent the values of y and Od during timestep J.
Volumetric microfracture density mrR3o
Now, the microfracture growth equation [11] can be used to derive expressions for the cumulative microfracture density distribution functions MFP3O and MFP32.
The volumetric microfracture density can be described by a cumulative population distribution function MFP3o(rN,tN), which gives the number of microfractures of radius GN or greater per m3 volume at the current time tN.
Various studies have shown that the initial microfracture population MFP3o(r,0) will follow a power law population distribution function:
[13] mRR3O(G, 0) = Br-c where B and c are constants; B is the initial microfracture density coefficient and c is the microfracture size distribution exponent; studies have shown that typically c=2. This can be combined this with [11] to derive an expression for the evolution of the volumetric microfracture density distribution function over time:
Figure imgf000028_0002
where GN represents a microfracture radius at time tN at the end of timestep N.
Note that the volumetric microfracture density given by [14] will include those microfractures pF that have reached the layer boundaries and become layer-bound macrofractures MF. To calculate the population of microfractures with r<h/2, the layer-bound macrofracture population must be subtracted from [14]
Mean linear microfracture density m R32
The total area of all microfractures in a unit volume, MFP32(0,t), can be calculated by integrating the number of microfractures of a given radius r times the area of a microfracture with radius r for radii between 0 and h/2 (the maximum radius of a microfracture in the layer). The total area of microfractures per unit volume is equivalent to the number of microfractures that will be encountered along a linear transect perpendicular to the fractures and, thus, can also be called the mean linear microfracture density.
The number of microfractures with radius between r and r+dr can be obtained by differentiating the cumulative volumetric microfracture density function MFP o(r,t) with respect to r. Therefore, the total area of all microfractures pF in a unit volume can be expressed as
Figure imgf000029_0001
Assuming b>2, dMFP o/dr at time tN can be obtained by differentiating [14] with respect to r:
<^Fp30 (r.tjy)
[16]
dr
Figure imgf000029_0002
So [15] becomes h s i N. -( C+1)
[17] mRR32 (0, tjv) = TTCB J0 2 gb+ 1 (r b + GN) dr
Unfortunately there is no analytical solution to [17]. However the limits of the integral (0 and h/2) are well defined and constant, so it is possible to obtain an approximate solution for [17] numerically, e.g.:
Figure imgf000030_0001
In most cases, a reasonably accurate result can be obtained with 1=10. The cumulative linear microfracture density as a function of microfracture radius, MFP32(rN,tN), can be calculated by adjusting the lower limit of the numerical integration in [18].
Modelling macrofractures In the following, it is examined how of circular microfractures pF can develop into layer-bound macrofractures MF, and expressions for the rate of macrofracture nucleation and growth are derived. These are used to calculate the cumulative density distribution functions describing the layer-bound macrofracture population, again ignoring fracture interaction.
Once the radius of a propagating microfracture pF reaches h/2, the fracture will span the entire layer and can no longer propagate vertically. It can, however, continue to propagate laterally, i.e. along the layer. It can thus be approximated as a long rectangular layer-bound fracture, which will be labelled as a macrofracture MF. Macrofractures MF will propagate independently at both lateral tips LT until propagation ceases for some reason (for example they intersect an orthogonal macrofracture MF or they enter the stress shadow of another parallel macrofracture MF). The key parameter characterising macrofracture size is thus the macrofracture length L.
Macrofracture propagation rate
It is assumed that the macrofracture length L is significantly greater than the layer thickness (and thus macrofracture height) h. The stress intensity Kf at the fracture tip LT is thus dependent only on h and not on L:
Figure imgf000031_0001
where the driving stress Od is as defined in [6]. [19] can be used to derive an equation for the lateral fracture propagation rate that is similar to [7]:
Figure imgf000031_0002
It should be noted that the factor of 2 is included in [20] to represent a macrofracture propagating at both tips, as each tip will propagate at half the rate given by [20] .
The key point to note about [20] is that it is independent of macrofracture length L. Therefore all macrofractures MF will propagate at the same rate, regardless of size, although this rate will vary with time as the driving stress Od changes. It is therefore not necessary to solve a differential equation to calculate macrofracture length L, which makes the calculation of the cumulative population distribution functions much easier. It is, however, necessary to consider that macrofractures MF will nucleate continuously through time, as microfractures mR grow to the maximum radius h/2.
In the following, it is shown how to calculate the cumulative population distribution functions when all macrofractures MF remain active and do not cease propagating for any reason. Further below, the method is extended to calculate active
(propagating) and static (non-propagating) macrofracture distribution functions when macrofractures MF can cease propagating (e.g. if they intersect with an orthogonal macrofracture MF or they enter the stress shadow of another parallel macrofracture MF).
Volumetric macrofracture density MFP30 The volumetric macrofracture density can be described by a cumulative population distribution function MFP3o(LN,tN), which gives the number of macrofractures of length LN or greater per m3 volume at the current time tN. MFP3o(l-N,tN) can be calculated by extrapolating back to the time tm at which a fracture with current length LN nucleated, at which time it had the length Lm=0. Because all macrofractures MF propagate at the same rate, the number of macrofractures of length LN or greater at time tN is equal to the total number of macrofractures at time tm, which in turn is equal to the total number of microfractures with radius h/2 or greater at time tm, MFP3o(h/2,tm):
Figure imgf000032_0001
If the driving stress Od varies through time then so will the macrofracture propagation rate dL/dt, and this must be taken into account when determining tm. If time is split into a series of timesteps J, each of duration tj-tj.i and with driving stress Odj, then the maximum macrofracture propagation distance 6LN,M between the end of timestep M and the end of timestep N will be given by
[23] SLn,M åj M+l[2aMFadJS (tj tj-l)] If tN lies at the end of timestep N and tm lies within timestep M (tM-i£tm<tM) then
Figure imgf000033_0001
It is helpful to define a function LN,M such that tj-l)]
Figure imgf000033_0003
Substituting [27] into [14] gives:
Figure imgf000033_0002
where timestep M is defined such that
[29] 2HMFAW,M £ LN < 2aMrALί Lί-1
Mean linear macrofracture density MFP32
The mean linear macrofracture density MFP32 can be calculated in the same way as for microfractures pF, by differentiating the volumetric macrofracture density MFP30 with respect to macrofracture length L to obtain the number of macrofractures MF with length between L and L+dL, multiplying the result by fracture area hl_, and integrating across the full range of possible macrofracture lengths L. For macrofractures MF, however, the fracture height is constant, so the macrofracture area is linearly proportional to the macrofracture length L. This simplifies the calculation of MFP32, as there is an analytical solution to the resulting differential equation:
Figure imgf000034_0001
[30] can be solved using the technique of integration by parts, to give
Figure imgf000034_0002
where timestep K represents the timestep of nucleation of a macrofracture with current length LN, defined such that
[32] 2aMFAW K· < LN < 2aMrLLί K·_1
[31] can be simplified by cancelling terms across the timesteps, to obtain
Figure imgf000035_0001
By setting LN=0, an expression for the total mean linear density of all macrofractures can also be derived:
Figure imgf000035_0002
Calculating macrofracture porosity and stress shadow volume
[33] can be converted into an expression for the total porosity <1>MF of a set of Mode 1 dilatant macrofractures 101 if the mean fracture aperture is known. It has been shown that if a cross section is taken of a laterally extensive elastic Mode 1 fracture 101 of height h (i.e. a 2-dimensional plane strain approximation), it will dilate with an elliptical profile, with aperture such that the total volume of a macrofracture of length L, VMF(L), will be given by
Figure imgf000035_0003
Therefore, the total macrofracture porosity at time tN, CPMF^N), will be given by
Figure imgf000035_0004
The total fracture porosity <1½F for a set of vertical Mode 1 dilatant macrofractures 101 also represents the total horizontal extensional strain accommodated by displacement on the fracture set, MFCh(t). Similarly, for a set of Mode 2 shear macrofractures 102, the total horizontal extensional strain accommodated by displacement on the fracture set is represented by the total volumetric heave in a unit volume PMR(ΪN). It has been shown that the displacement of a Mode 2 fracture 102 follows an identical profile to that of a Mode 1 fracture 101 , and hence the total volumetric heave can be calculated by simply adding an additional term to the formula for total macrofracture porosity [36], to account for the fracture dip w:
Figure imgf000036_0001
Therefore, an expression can be derived for the total horizontal strain
accommodated by a set of macrofractures MF in terms of the mean linear macrofracture density MFP32:
Figure imgf000036_0003
where the minus signs indicate extensional strain, and w is the fracture dip.
Fractures are often surrounded by a stress shadow, which comprises the zone in which a change in the applied strain is accommodated by an elastic displacement on the fracture, rather than by a change in the elastic strain (and hence stress) in the rock mass surrounding the fracture. The total width of the stress shadow is given by the change in fracture displacement divided by the change in bulk rock elastic strain. Therefore, [38] can be combined with an expression for the horizontal stress (i.e. the driving stress) as a function of applied horizontal strain, to obtain an expression for the stress shadow width around a macrofracture. The relationship between driving stress and applied horizontal strain in an elastic rock mass can be obtained by combining [6] with the 3-dimensional form of Hooke’s Law, giving
Figure imgf000036_0002
for a Mode 1 macrofracture 101 perpendicular to the minimum horizontal applied strain Shmin, and
[40] sά = (sin
-(sin w c
Figure imgf000037_0001
for a Mode 2 macrofracture 102 striking perpendicular to Chmin, with Chmax representing the maximum horizontal strain, ov’ the effective vertical stress, and E and v the Young’s Modulus and Poisson’s ratio of the bulk rock mass.
The mean stress shadow width W around individual macrofractures is thus given by
Figure imgf000037_0002
for Mode 1 macrofractures 101 with mean aperture a, and
Figure imgf000037_0003
for Mode 2 shear macrofractures 102 with mean displacement d.
[38] can also be combined with [41] and [42] to obtain expressions for the total macrofracture stress shadow volume IJJMF as a proportion of total rock volume:
Figure imgf000037_0004
for Mode 1 dilatant fractures 101 , and
[44]
(sin w cos w + m G sin2 w)MRR32 (0, ΐL,) for Mode 2 shear fractures 102. Note that, since stress shadows cannot overlap, the total stress shadow volume is simply the sum of the volume of the stress shadows around every macrofracture.
Active and static fractures
So far, it has been assumed that microfractures will all continue propagating until they reach radius r=h/2, and that macrofractures will continue propagating indefinitely, as long as there is a driving stress. However in practice fractures may cease propagating at any time, if they intersect an orthogonal fracture, or if they fall into the stress shadow of another parallel fracture. The equations for the cumulative fracture population distributions P30 and P32 must therefore be modified to take account of this.
This can be done by dividing the microfractures and macrofractures at any time t into two populations: active fractures (apF and aMF), which are still propagating, and static fractures (spF and SMF), which have ceased propagating. All fractures are initially considered to be active. However, active fractures can cease propagating and become static fractures at any time. Static fractures cannot revert to active fractures, and will therefore remain the same size indefinitely. It is important however to note that static fractures can still accumulate displacement and accommodate applied strain. By considering the different mechanisms by which microfractures and
macrofractures may cease propagating, the probability that an active fracture will be deactivated and become a static fracture can be calculated. This may take one of two forms: the instantaneous probability of fracture deactivation represents the probability that a specific fracture will cease propagating within a very short time increment dt, while the cumulative fracture activation probability represents the probability that an initially active fracture will still be active at the end of a longer time interval, e.g. a timestep. If fracture deactivation is a random process, so that the instantaneous probability of fracture deactivation q is independent of time, then the cumulative fracture activation probability Q will follow an exponential decay curve of the form
[45] 0 (t) = e-
However, in some cases fracture deactivation will be a regular or semi-regular process, where the instantaneous probability of fracture deactivation increases with time. In these cases, it is easier to calculate the probability that a fracture will be deactivated over a longer time interval (i.e. the cumulative activation probability), and then divide this by the time interval to obtain a pseudo-instantaneous or mean probability of fracture deactivation.
The volumetric density of active fractures (VP30 and 8 MFP3O) can be calculated by combining the cumulative fracture activation probability with the expressions for MFP 30 or MFP30 derived above. If the resulting active fracture volumetric densities, VP30 and 8 MFP3O, are then multiplied with the instantaneous fracture deactivation probability (or the mean fracture deactivation probability), the rate at which active fractures cease propagating and become static fractures is obtained. Integrating these deactivation rates through time will therefore give expressions for the volumetric density of static fractures, VP30 and S MFP3O. Finally, the previous procedure of differentiating the P30 functions can be repeated to determine the number of fractures with a given area, multiplying this by the fracture area, and integrating the resultant functions across the range of fracture sizes to obtain functions for the volumetric active and static microfracture ratios and mean linear active and static macrofracture densities, VP32, 8MFP32, VP32, and SMFP32.
One advantage of this procedure is that the calculation of the fracture deactivation probabilities can be separated from the derivation of the active and static fracture population functions VP30, 3MFP3O, VP30, SMFP3O, VP33, 8MFP33, VP32, and P32: the fracture population functions include expressions for the fracture deactivation probabilities but are independent of how these are calculated. This is important because, since fractures cease propagating as a result of interaction with other fractures, the fracture deactivation probabilities will change through time as the fracture populations grow. By calculating the fracture deactivation probabilities independently from the fracture population functions, the timestep approach can be used to include this feedback: the fracture deactivation probabilities at the start of each timestep can be calculated, based on the fracture populations at the end of the previous timestep, and then these can be used to calculate the evolution of the fracture population distribution functions during the subsequent timestep.
This approach, however, is not valid if the average lifetime (or half-life) of an active fracture is shorter than the timestep duration. Unfortunately, this is likely to be the case in the later stages of the growth of the fracture network, when the probabilities of fracture interactions are high. However, in this situation equations describing the “residual” fracture populations present in these timesteps can be derived as an equilibrium population, where individual fractures are short-lived but the rate of fracture nucleation is equal to the rate of fracture deactivation.
Calculating the fracture deactivation probabilities
Three mechanisms by which fractures may cease propagating will be considered, one of which applies to microfractures and two of which apply to macrofractures:
• Microfractures will cease propagating if they fall into the stress shadow of a nearby parallel propagating macrofracture.
• Macrofractures will also cease propagating if their stress shadow interacts with the stress shadow of another parallel macrofracture.
• Macrofractures will cease propagating if they intersect an orthogonal or oblique macrofracture.
Probability of microfracture deactivation due to interaction with macrofracture stress shadows
To simplify the calculation of the microfracture deactivation probability, some assumptions about microfracture behaviour are made: • The microfractures are small compared to the macrofractures, and propagate at a slower rate.
• Therefore a microfracture will cease propagating if the centre of the
microfracture falls within the stress shadow zone of a parallel macrofracture. · The microfracture will have no effect on the macrofracture, which will
continue propagating as before.
• The microfractures have insufficient size or density to develop a significant volume of stress shadow themselves. Interaction with other microfractures is therefore ignored.
If it is assumed that all microfractures are active at time t=0, and that they are evenly distributed throughout the layer, then the proportion of microfractures that will have been deactivated at any later time tN will be given by the ratio of the macrofracture stress shadow volume to the total volume of the layer at time tN. Therefore, the proportion of microfractures still active at time tN, i.e. the cumulative microfracture activation probability QN, is given by:
[46] QN— 1 MF ( N) The probability that a microfracture that is not currently within a macrofracture stress shadow will fall into a macrofracture stress shadow within a very short time interval dt, i.e. the instantaneous probability of microfracture deactivation, q, is purely a function of the rate of increase of the total macrofracture stress shadow volume IJJMF:
Figure imgf000041_0001
Note that combining [46] and [47] results in a differential equation for QN that is consistent with [45]. It can also be seen from [47] that the instantaneous probability of microfracture deactivation will increase continuously, as the total macrofracture stress shadow volume IJJMF increases. However a mean rate of microfracture deactivation qwi for a given timestep M can also be defined by dividing the total probability of microfracture deactivation during timestep M by the timestep duration:
Figure imgf000042_0001
One final issue is that [46] assumes that the microfractures themselves have no stress shadow. Therefore if the probability that a microfracture will grow into macrofracture, with its own stress shadow is to be calculated, [46] must be modified to define an exclusion zone around each macrofracture, equal to twice the stress shadow width W, so that if a macrofracture nucleated within the exclusion zone the two stress shadows will overlap. Since the exclusion zones may potentially overlap, the volume of the exclusion zones surrounding the macrofractures,
Figure imgf000042_0002
must be calculated as described below. This results in an effective cumulative microfracture activation probability Q'N:
Figure imgf000042_0003
Probability of macrofracture deactivation due to stress shadow interaction
Like microfractures, layer-bound macrofractures may cease propagating due to stress shadow interaction. However there are several key differences:
• In order for the two macrofracture stress shadows to interact, at least one of the macrofracture tips must be active and propagating towards the other.
• All macrofractures have stress shadows, and stress shadow interaction will occur when the stress shadows around two parallel macrofractures overlap.
• When this happens, both of the macrofracture tips will be deactivated.
The process of macrofracture stress shadow interaction is illustrated schematically in Fig. 3. An exclusion zone can be defined around each macrofracture, equal to twice the stress shadow width W, so that if another macrofracture tip enters that exclusion zone, the two stress shadows will overlap and both fracture tips will be deactivated. Unlike the stress shadows, however, the exclusion zones around the macrofractures can partially overlap. If the macrofractures are randomly distributed, but in such a way that that their stress shadows do not overlap, it can be shown that the total exclusion zone volume XMF will be given by
Figure imgf000043_0001
If a very short time interval dt is considered, so that the increase in the total macrofracture exclusion zone volume during dt is small, then a propagating macrofracture tip LTA of a first macrofracture MFA with a stress shadow 201 will interact with the stress shadow 301 around a second parallel macrofracture MFB with exclusion zone 302 if the macrofracture tip LTA lies within the so-called stress shadow interaction box 303 projected ahead a macrofracture tip LTB of
macrofracture MFB, as illustrated in Fig. 3. This stress shadow interaction box 303 will have width 2W, where W is the mean stress shadow width, and length LSSIB equal to the macrofracture tip propagation rate (i.e.
Figure imgf000043_0002
Odb, cf. [20]) times dt, if the macrofracture tip LTB is static, or twice that if both the macrofracture tips LTA, LTB are propagating towards each other. The probability of stress shadow interaction between the two macrofractures MFA, MFB within time interval dt is therefore given by the volume of the stress shadow interaction box 303 as a proportion of the total “clear zone” volume (i.e. the volume not lying within a macrofracture exclusion zone, given by 1 -XMF). Furthermore, the probability that a propagating macrofracture tip LT will experience stress shadow interaction with any other macrofracture MF within time interval dt will be given by the total volume of the stress shadow interaction boxes 303 of all parallel macrofracture tips LT facing in the opposite direction, whether static or dynamic. This will be proportional to the volumetric macrofracture density MFP3O.
Since macrofractures MF become deactivated due to either stress shadow interaction or intersection, all static macrofracture tips LT must lie either directly adjacent to the stress shadow zone of another parallel macrofracture MF, or else on an orthogonal macrofracture MF as is the case for LTB in Fig. 4. In this case, in which the second macrofracture MFB has been deactivated due to intersection with a third orthogonal macrofracture MFc, the propagating macrofracture tip LTA will be deactivated due to intersection with the orthogonal macrofracture MFc before it reaches the exclusion zone 302 of the static macrofracture tip LTB. Thus, such “obstruction effects” will reduce the probability of stress shadow interaction between the propagating and static half-macrofractures, and reduce the effective size of the stress shadow interaction boxes.
Therefore, a function x(ΪN) is defined, representing the total rate of increase of the combined length of all stress shadow interaction boxes facing in one direction at time tN, weighted to take account of obstruction effects:
[51 ] x(ίN)
Figure imgf000044_0001
where SM MFP3O and SU MFP3O represent the volumetric density of macrofractures that are deactivated due to stress shadow interaction and intersection respectively, and uV and uV are the modifiers to take account of obstruction effects at the static macrofracture tips. In general, the mean overlap of exclusion zones of two macrofractures terminating due to stress shadow interaction will be 75%, so it can be assumed that uV = 0.25. In the case of fracture intersection, unless the intersection point lies close to the tip LTc of the orthogonal macrofracture MFc, the propagating macrofracture tip LTA will always intersect the orthogonal macrofracture MFc before it reaches the stress shadow 301 of the parallel macrofracture MFB. Therefore, it can approximated that uV « 0 (although it is possible to calculate a more precise value if the macrofractures MF are short).
Then two new functions can be defined:
Figure imgf000044_0002
The total volume of stress shadow interaction boxes for the interval dt will therefore be given by dxp/dt dt. If the timestep N is reasonably short, it can be assumed that HJMF, XMF and x will remain constant throughout the timestep and the values calculated at the end of the previous timestep can be used, so that YMR,N=YMR(IN-I), XMF,N=XMF(tN-i) and x N= x (ΪN ) The instantaneous probability of macrofracture stress shadow interaction FI,I,N will then be given by
Figure imgf000045_0001
Since macrofracture stress shadow interaction is a random process (i.e. FI,I,N is independent of the time that the macrofracture MF has been active), the cumulative macrofracture activation probability for stress shadow interaction FI,I,N is given by an exponential function:
Figure imgf000045_0002
Probability of macrofracture deactivation due to intersection
As mentioned above, layer-bound macrofractures MF will also cease propagating if they intersect a perpendicular-striking macrofracture MF, as long as the frictional traction on the surface of the intersected fracture is sufficiently low to allow slip. It is assumed that this is always the case, so the probability of macrofracture intersection is determined only by the propagation rate of the intersecting macrofracture tip LT and the population of perpendicular-striking macrofractures MF.
However, the probability of macrofracture intersection is dependent not only on the density of the perpendicular-striking macrofractures MF, but also on their spatial distribution. If the perpendicular macrofractures MF are randomly distributed, fracture intersection will be a random process, governed by an exponential probability function. However, if the perpendicular macrofractures MF are evenly spaced (for example if their distribution is controlled by stress shadows), then fracture intersection will be a regular or semi-regular process. In reality, many fracture sets will be somewhere in between these two end-members: There will be a minimum separation between the macrofractures MF, determined by the stress shadow width, but there will also be larger gaps that are randomly distributed.
Macrofracture intersection then becomes a semi-regular process that is especially hard to model, as it will include elements of both regular and random behaviour.
A general expression for the probability of macrofracture intersection that is valid for any distribution of macrofractures MF can be derived by considering a macrofracture tip LT from a first set of macrofractures MF propagating towards a second set of perpendicular-striking macrofractures MF. The macrofractures MF of this second set are surrounded by stress shadows of width W, which cannot overlap. W may take any value from 0 (for a fully randomly distributed fracture set) to 1/J MFP32 (for an evenly spaced set), where J MFP32 is the mean linear density of the second set of macrofractures.
If a short timestep of duration At is taken, the propagating fracture tip LT will intersect a macrofracture MF of the second set during At if that macrofracture MF lies within the distance DI that the macrofracture tip LT of the first set will propagate in the time period At. A set of“intersection boxes” with width DI can therefore be defined, extending towards the propagating macrofracture tip LT from each of the perpendicular-striking fractures MF. The probability of fracture intersection during the timestep is thus equal to the total volume of all intersection boxes for the second set of macrofractures MF as a proportion of the total rock volume.
The total volume of an intersection box attached to a specific macrofracture MF of the second set with length L is given by LhAI. If the width DI of the intersection boxes is less than the stress shadow width W then none of the intersection boxes can overlap, and the total volume of all the intersection boxes is given by J MFP AI.
However, if the width of the intersection boxes is greater than the stress shadow width W, those parts of the intersection box that lie further from the macrofracture MF than the stress shadow width W may overlap the intersection boxes around other macrofractures MF. In this case, a“one-sided proximity zone” formula must be used to calculate the total intersection box volume, similar to that used to calculate total exclusion zone volume in the previous section.
Combining these two forms gives a general expression for the probability of macrofracture intersection during a timestep N with duration DΪN, i.e. the cumulative macrofracture activation probability for fracture intersection <PI,J,N:
Figure imgf000047_0002
Hill is the integer part of i, rounding down, Ϊ is the fractional part (i =i-llill), and
Figure imgf000047_0001
where 'odN is the driving stress acting on the first set of macrofractures. Note that if W=0 (i.e. there are no stress shadows) then i®¥ and [58] can be reduced to
Figure imgf000048_0001
Finally, if the probability that an initially active macrofracture will be deactivated during timestep N is divided by the timestep duration, the mean probability of macrofracture intersection F N can be calculated:
Figure imgf000048_0002
Combining the macrofracture deactivation probabilities
Since the two macrofracture deactivation mechanisms are independent, the probability FN that a macrofracture which is active at the start of a timestep N will remain active at the end of the timestep can be found by multiplying the cumulative macrofracture activation functions for stress shadow interaction and for intersection:
[63] OJV = F/,/,LίF/jw
The instantaneous probability of macrofracture deactivation for the timestep FN is much harder to calculate since, as macrofracture deactivation can occur by a combination of random and regular processes, it will vary throughout the timestep in a complex manner. However, the mean probability FN of macrofracture deactivation can easily be calculated by dividing the probability of macrofracture deactivation during the timestep by the timestep duration:
Figure imgf000048_0003
As it will be seen in the following sections, it is also useful to define the cumulative probability FN,M that a macrofracture that is active at the end of a timestep M will remain active at the end of a later timestep N:
[65] F«,M — FM+1 X FM+2 X " X F JV
Calculating active and static microfracture populations
The microfracture deactivation probabilities can be combined with the cumulative microfracture population functions derived above to obtain expressions for the active and static microfracture populations.
Volumetric microfracture density
The volumetric active and static microfracture density functions VP30 and VP30 can be calculated by modifying the equation [14] for volumetric microfracture density.
The active microfracture density VP30 is obtained by multiplying the total microfracture density by the cumulative microfracture activation probability:
Figure imgf000049_0001
To calculate the volumetric static microfracture density S MFP3O, an expression for the rate of microfracture deactivation must be derived and integrated through time from t=0 to t=tN. The rate of microfracture deactivation is equal to the active microfracture population times the instantaneous microfracture deactivation rate. However, as the cumulative microfracture activation probability can only be calculated at the end of each timestep, a constant value for Q must be used throughout each timestep, reflecting the macrofracture population at the end of the previous timestep.
Therefore, the integral must be discretised into a sum of terms representing individual timesteps. Also, the mean probability of microfracture deactivation q« rather than the instantaneous probability of microfracture deactivation q« should be used when calculating the rate of microfracture deactivation during a timestep K. This gives
Figure imgf000050_0001
Summing the integral of [67] across all timesteps up to and including timestep N gives
[68]
Figure imgf000050_0002
Note that [67] and [68] include microfractures that have reached the layer boundaries and become layer-bound macrofractures. To calculate the population of microfractures with r<h/2, the layer-bound macrofracture population must be subtracted from [67], and an extra term must also be included in [68] to account for microfractures that are deactivated as they reach r=h/2 because they lie in a macrofracture exclusion zone:
Figure imgf000050_0003
Mean linear microfracture density
The mean linear active and static microfracture densities VP32 and VP32 can be calculated using the same method as previously: by differentiating the equations for volumetric microfracture density to determine the number of microfractures with a given radius r, multiplying this by the microfracture area or volume, and integrating the resultant function from r=0 to r=h/2.
As noted above, there is no analytical solution to the integral equations for the mean linear microfracture density MFP32 when b¹2, but a reasonable approximation can be obtained by numerical integration:
Figure imgf000051_0001
Calculating the linear static microfracture density VP32 is more complex as the summations must be differentiated across timesteps given in [68] with respect to r, multiplied by fracture area, and then reintegrated numerically across the range of microfracture radii from r=0 to r=h/2, to get:
Figure imgf000051_0002
Calculating active and static macrofracture populations
Expressions for the active and static macrofracture populations can be obtained in a similar manner: combining the macrofracture deactivation probabilities with the cumulative microfracture population functions derived above.
However when modelling macrofractures MF, it must recognised that the two macrofracture tips LT will propagate independently, and each can become deactivated independently of the other. The easiest way to account for this is to model half-macrofractures ME rather than macrofractures MF. A half-macrofracture ME represents the portion of a macrofracture MF between the nucleation point and one lateral tip LT; each macrofracture MF therefore comprises two conjoined half- macrofractures ME- Therefore, the cumulative macrofracture population functions derived above can easily be converted into cumulative half-macrofracture population functions, expressed in terms of half-macrofracture length I, by replacing the initial microfracture density coefficient B with 2B and replacing the macrofracture propagation rate dL/dt, given by [20], with the half-macrofracture propagation rate dl/dt, given by
Figure imgf000052_0001
Volumetric half-macrofracture density
The volumetric fracture density functions, a MFP3o(l,t) and s MFP3o(l,t), for active and static half-macrofractures ME can be calculated by modifying [28]. However, since these calculations are complex, it may be helpful to represent the growth of the fractures graphically. Fig. 5 illustrates the growth through time of half-macrofractures ME nucleating at different times, specifically at the end of each timestep (the time tx indicating the end of timestep X, timestep N being the current timestep) as a series of lines on a graph of half-macrofracture length against time. The solid lines represent the growth of half-macrofractures ME that nucleate at the timestep boundaries and remain active until the end of timestep N. Of course, macrofractures MF may nucleate at any time, and the dashed line indicates the growth trajectory of a half-macrofracture ME that nucleated in the middle of timestep 7. Note that within each timestep, all half-macrofractures ME grow at the same rate, but that this rate can vary between timesteps (in the example shown in Fig. 5, the rate decreases through time, but this does not have to be the case). Thus, at any time t, the time of nucleation of a half-macrofracture ME of length I can be determined by tracing back along the growth curve.
Active volumetric half-macrofracture density 3 MFP3O
Calculating the volumetric density of active half-macrofractures ME is more complicated than calculating the volumetric density of active microfractures pF, since half-macrofractures ME of different lengths I will have nucleated at different times. The probability that a half-macrofracture ME is still active (expressed in terms of the cumulative macrofracture activation probability F) will depend on its nucleation time tm, and therefore also on its length I. The probability that the seed microfracture pF will have been deactivated before it could become a macrofracture MF must also be taken into account; this also depends on the time of macrofracture nucleation.
Therefore, [28] cannot simply be multiplied by a single cumulative fracture activation probability to obtain an expression for the cumulative volumetric density of active half-macrofractures ME. as is the case for microfractures pF. Instead, an expression for the number of currently active half-macrofractures ME as a function of half- macrofracture length I must be derived, and then integrated across the range of possible half-macrofracture lengths I.
To do this, [28] must first be modified to apply to half-macrofractures ME rather than to macrofractures MF, by substituting I for 2L and 2B for B. Then, the maximum density of half-macrofractures ME with length between I and l+dl, if there is no fracture deactivation, is calculated by differentiating the modified version of [8] with respect to I. The nucleation time tm for a half-macrofracture ME of current length I must also be determined by extrapolating back along the half-macrofracture growth trajectory (the solid lines in Fig. 5). Then, the maximum number of half- macrofractures ME can be multiplied by the cumulative fracture activation probability for a fracture nucleating at time tm, and by the effective microfracture activation probability for the time tm, to obtain the density of active half-macrofractures ME with length between I and l+dl. Finally, this function must be integrated between the specified lower half-macrofracture length limit IN and the current maximum half- macrofracture length, i.e. the length of a macrofracture that nucleated at time t=0, which is given by OMFAN.O. However, this calculation can be simplified by discretising the cumulative activation probabilities to the timestep boundaries. Hence FN,M-I is defined as the probability that a half-macrofracture ME that nucleated at the start of timestep M is still active at the end of timestep N, and QM-I is defined as the probability that a given
microfracture pF is active at the start of timestep M. These discretised cumulative activation probabilities are then used to calculate the active densities of all half- macrofractures ME that nucleated at any time during timestep M. When this assumption is made, the integral equation described above can be replaced with a sum across timesteps, so that [28] becomes
Figure imgf000054_0001
where timestep J is defined such that
Figure imgf000054_0002
This is illustrated graphically in Fig. 5. To calculate the cumulative volumetric density of all active half-macrofractures ME with length greater than IN at time tN using continuous fracture activation probabilities, the current active density of all half- macrofractures ME with lengths greater than IN, would be integrated along the vertical arrow at the right side of the diagram. If, however, the fracture activation probabilities are discretised by timestep, this integration is converted into a sum of discrete terms representing each timestep from timestep 1 , when the longest currently active half-macrofractures ME nucleated, to timestep J, when half- macrofractures ME with length IN nucleated. The boundaries of these terms are the points at which the solid half-macrofracture growth trajectories cross the integration arrow; between these boundaries, the cumulative fracture activation probabilities are uniform. A semi-continuous form of [73] for any time tn within timestep N can also be derived:
Figure imgf000055_0001
This will be used when calculating the static volumetric half-macrofracture density. Note that although the term for maximum fracture density is continuous with respect to time in [75], the cumulative activation probabilities remain discretised.
Static volumetric half-macrofracture density S MFP3O
The cumulative distribution functions for the volumetric density of static half- macrofractures ME at time tN can be calculated in the same way as was done for the volumetric density of static microfractures pF: the volumetric density of active half- macrofractures ME (calculated using [75]) is multiplied with the mean probability of half-macrofracture deactivation FK, to obtain an expression for the rate of deactivation of half-macrofractures ME of length IN or greater, dS MFP3o(lN,t)/dt. Then, this is integrated through time from t=0 to t=tN to obtain S MFP3o(lN,tN). This is illustrated in Fig. 6. To calculate the static half-macrofracture density, the active half-macrofracture density times the half-macrofracture deactivation rate is integrated with respect to time, along the horizontal arrow in Fig. 6. However, the active half-macrofracture density is itself an integral with respect to half-macrofracture length I, so overall it is necessary to integrate across the entire hatched area in Fig. 6. J represents the timestep in which an active half-macrofracture ME must have nucleated in order to have length IN at a specific time t. Moving along the integration path, the value of J will change whenever a half-macrofracture growth trajectory line is crossed. This will generally not coincide with the timestep boundaries. In Fig. 6, for instance, timestep K=9 has multiple J values, whereas in timestep K=1 1 , the J value is constant.
As with microfractures pF, this calculation can be simplified by discretising it into a series of timesteps K, each with a constant mean probability of half-macrofracture deactivation FK. This gives
Figure imgf000056_0001
Now, the semi-continuous form that was derived for the volumetric density of active half-macrofractures [75] can be used. Inserting [75] into [76] and integrating across each timestep K gives
Figure imgf000056_0002
where timestep J is such that
Figure imgf000056_0003
However, it must be remembered that the active half-macrofracture volumetric density equation [75] is itself an integral with respect to half-macrofracture length I, albeit discretised to a sum of terms representing the half-macrofracture nucleation timesteps. The last term in this sum, which represents the nucleation timestep of a half-macrofracture of length IN, varies continuously with time. This is illustrated in Fig. 6, where the horizontal arrow represents the path of integration through time, the vertical thin lines represent the boundaries of the discretised values of the mean fracture deactivation probability FK (dependent on the timestep of fracture deactivation K in [77]), and the vertical arrows represent the integration along these boundaries and along the boundaries of the nucleation timestep of a half- macrofracture of length IN (i.e. the timestep J in [75]). These arrows are made thicker for the timesteps referred to below (K=9 and K=1 1 ). As it can be seen in Fig. 6, the two sets of boundaries do not generally coincide, so there may be multiple timesteps J that satisfy [78] for different values of t between -i and (e.g. timestep K=9 in Fig. 6). This makes the definite integral across timestep K in [77]
complicated.
The easiest way to solve [77] is to integrate each of the individual timesteps M in the inner sum (including the final term for M=J) separately across timestep K. Thus, [77] can be expressed as a second inner sum over timesteps from M=1 to M=K, where M represents the timestep of nucleation and K represents the timestep of deactivation of a subset of static half-macrofractures of half-length IN or greater:
Figure imgf000057_0001
where <8 MER3O>M represents the terms for individual timesteps M in the sum in [77] (including the final term for M=J), and <|8 MRR3O>K,M represents the integral of these terms across timestep K. The form of <8 MFP3O>M will vary depending on whether M<J (when it will be independent of time t), M=J (when it will be time-dependent), or M>J (when it will be 0). The <|8 MRR3O>K,M terms in the inner sum can therefore be one of six different types, depending on whether the M=J for some, all or none of timestep
K:
Type 1 terms
For timesteps M where IN £ OMFAK-I,M and IN < OMFAK,M (e.g. timesteps M=1 , M=2 and M=3 for K=9 in Fig. 6) then M<J for all values of t between k-i and k, so only the time independent component of the integral is relevant. These terms are therefore given by
Figure imgf000058_0001
Type 2 terms
For timesteps M where 0MFAK-I,M < IN £ 0MFAK-I,M-I and IN < 0MFAK,M (e.g. timestep M=4 for K=9 in Fig. 6) then M=J for values of t between k-i and
tK-(aMFAK,M-lN)/(aMFadKb), but M<J for values of t between ΐk-(aMrLk,M-ΐN)/( aMFadK b) and k. Therefore, both the time dependent component of the integral (when M=J) and the time independent component of the integral (when M<J) must be included to obtain the term:
Figure imgf000058_0002
Type 3 terms
For timesteps M where OMFAK-I,M-I < IN < OMFAK,M (e.g. timestep M=5 for K=9 in Fig.
6) then M>J for values of t between -i and tK-(aMFAK,M-i-lN)/( aMFadKb), M=J for values of t between tK-(aMFAK,M-i-lN)/( aMFadK b) and tK-(aMFAK,M-lN)/( aMFadK b), and M<J for values of t between tK-(OMFAK,M-lN)/( OMFadKb) and . The interval where M>J can be ignored, but both the time dependent component of the integral (when M=J) and the time independent component of the integral (when M<J) must be included to obtain the term:
Figure imgf000059_0001
Type 4 terms
For timesteps M where 0MFAK,M £ IN < 0MFAK,M-I and IN > 0MFAK-I,M-I (e.g. timestep M=6 for K=9 in Fig. 6) then M>J for values of t between -i and
tK-(aMFAK,M-i-lN)/(aMFadKb) and M=J for values of t between tK-(aMFAK,M-i-lN)/( aMFadK b) and . The interval where M>J can be ignored, but the time dependent component of the integral (when M=J) must be included to obtain the term:
Figure imgf000059_0002
Type 5 terms
It is also possible to have a timestep M where 0MFAK-I,M < IN £ 0MFAK-I,M-I but
IN ³ QMFAK,M (e.g. timestep M=7 for K=11 in Fig. 6). In this case, M=J for all values of t between k-i and . Therefore only the time dependent component of the integral needs to be included to obtain the term:
Figure imgf000060_0001
Type 6 terms
Finally, for timesteps M where IN ³ OMFAK,M-I (e.g. timestep M=7 for K=9 in Fig. 6) then M>J for all values of t between k-i and k. Therefore
Figure imgf000060_0002
and these timesteps cam be ignored.
An expression for the volumetric density of static half-macrofractures of length IN or greater at time tN can therefore be obtained by inserting the appropriate terms 1 to 6 into the double sum in [79]. Note that when IN = 0, a Type 1 term can be used for all timesteps M<K and a Type 5 term can be used for timestep M=K. The Type 5 term can be further simplified since IN=AK,M=AK-I,M-I =0. The total volumetric density of all static half-macrofractures can thus be written as
Figure imgf000060_0003
Mean linear half-macrofracture density
The mean linear density functions, a MFP32(l,t) and s MFP32(l,t), for active and static half- macrofractures ME can be calculated in the same way as MFP32(L,t): by
differentiating the appropriate P30 functions with respect to I to determine the number of half-macrofractures JV with half-length between I and l+dl, multiplying by Ih to obtain the total area of half-macrofractures ME with half-length between I and l+dl, and then integrating this expression from IN to the maximum possible half- macrofracture half-length ax.
Active mean linear half-macrofracture density 3 MFP32 The differential of MFP3o(l,t) in [30] can be modified to include the cumulative activation functions FN,M-I and O’M-I , and also to represent half-macrofractures ME instead of macrofractures MF. This gives
[87]
Figure imgf000061_0001
for a half-macrofracture ME that nucleated in timestep M. Since the cumulative activation probabilities FN,M-I and O’M-I are discretised to the timestep boundaries, they can be treated as constants when integrating across half-macrofractures ME that nucleated in the same timestep M. Since the full form of equation [31] derived above is already discretised by the macrofracture nucleation timestep, it can be modified relatively easily to obtain an expression for the mean linear density of active half-macrofractures ME of length IN or greater at time ΪN:
Figure imgf000061_0002
Figure imgf000062_0001
where timestep K represents the timestep of nucleation of a half-macrofracture with current length IN, defined such that
Figure imgf000062_0002
This time, however, [88] cannot be simplified by cancelling terms across the timesteps, since the cumulative activation probabilities will be different. An expression for the total mean linear density of all active half-macrofractures ME can be derived by setting IN=0:
Figure imgf000062_0003
Static mean linear half-macrofracture density S MFP32
The mean linear density for static half-macrofractures ME can be calculated in a similar way: by differentiating the volumetric density functions for static half- macrofractures ME with respect to length I to obtain an expression for the number of half-macrofractures ME with length between I and l+dl, multiplying this by half-length I and height h, and then integrating the resulting functions across a range of half- lengths, up to the maximum possible half-length ax. This is illustrated in Fig. 7.
In this diagram, the term for a specific combination of nucleation timestep M (in this case M=5) and deactivation timestep K (two examples are shown, K=9 and K=11 ), is integrated with respect to half-macrofracture length I (i.e. in the vertical direction). The nucleation timestep J of an active half-macrofracture ME with the exact length I at time t will vary with both t and I. Therefore, as I increases, the <.[3MFP3O>K,M term, representing the density of active half-macrofractures nucleated in timestep M, will vary between the six types described above, depending on whether M=J for some, all or none of timestep K. For any timestep K, the total value of the integral across all lengths between IN and Lax, and from -i to k, can be calculated by summing the integrals for each of the six types of the <.[3MFP3O>K,M term, over the intervals in which they are valid (note that some types may not be valid over any interval, for any given combination of K and M).
The functions for S MFP3O derived above are broken down into multiple terms based on the timestep of half-macrofracture nucleation M and the timestep of deactivation K. Each of these terms can be differentiated and reintegrated independently. Thus the differential of the static half-macrofracture volumetric density with respect to length is
Figure imgf000063_0001
and the mean linear density of static half-macrofractures s MFP32(lN, ) is given by
Figure imgf000063_0002
However as [92] is integrated across the range of possible half-macrofracture lengths from IN to Lax, the nucleation timestep for an active half-macrofracture ME with the exact length I (timestep J) will vary, and hence the form of the <[3MFP3O>K,M term will also vary between the six types described above, depending on whether M=J for some, all or none of timestep K. To account for this, the integral in [92] can be broken down into six components, each representing the range of half- macrofracture lengths I covered by a type T term <[3MFP3O>T, as defined above. This results in an extra summation:
Figure imgf000064_0001
where Linco and ax(T) represent the limits of the range of half-macrofracture lengths I covered by a type T term <.[3MFP3O>T, or IN if this is greater. This is illustrated for a specific timestep M=5 in Fig. 7. The integral components for each of the six types of terms <|3MFP32>T can be calculated as follows:
I ntegral component for Type 1 range
When I < OMFAK-I,M (and I < OMFAK,M) then [80] gives us
Figure imgf000064_0002
This component can therefore be ignored.
I ntegral component for Type 2 range
When QMFAK-I,M < I £ QMFAK-I,M-I and I < QMFAK,M then [81 ] gives us
Figure imgf000064_0003
Figure imgf000065_0001
Integral component for Type 3 range
When QMFAK-I,M-I < I < QMFAK,M (i.e. QMFAK,M > QMFAK-I,M-I) then [82] gives us
[97] (f MI¾q k,M
Figure imgf000066_0001
Figure imgf000067_0001
I ntegral component for Type 6 range
When I > QMFAK,M-I then [85] gives us
Figure imgf000068_0001
This component can therefore be ignored. The full integral from IN to lmax, for half-macrofractures nucleated in timestep M, is given by summing the integral components for each of the types. Note however that only those components where M=J for a half-macrofracture of length IN during some or all of timestep K (i.e. the components for types 2, 3, 4 and 5) will contribute to the integral. Combining these four term types for a specific timestep combination M,K gives us
Figure imgf000068_0002
[104] can be simplified if the value of IN relative to LK,M-I , LK,M, AK-I,M-I and AK-I,M is known. Again, six different types of terms can be defined, in this case reflecting the value of IN relative to the timesteps M and K (i.e. the lowest term component in the summation, represented by the horizontal arrow in Fig. 7). Type 1 terms
If IN £ QMFAK-I,M and IN < QMFAK,M (i.e. QMFAK,M > 0) then [104] can be reduced to
Figure imgf000069_0001
Type 3 terms
If QMFAK-I,M-I < IN < QMFAK,M (i.e. QMFAK,M > QMFAK-I,M-I) then [104] can be reduced to
Figure imgf000070_0001
Type 4 terms
If IN > QMFAK-I,M-I and QMFAK,M £ IN < QMFAK,M-I then [104] can be reduced to
[108] Type T [ / MF¾2 T]
Figure imgf000070_0002
Figure imgf000071_0001
Type 6 terms
If IN ³ QMFAK,M-I then [104] can be reduced to
Figure imgf000071_0002
The mean linear density of static half-macrofracture s MFP32(lN,tN) for any half-length IN can thus be obtained by inserting [105], [106], [107], [108], [109] or [1 10], as appropriate, into [93]. If IN = 0 then it can be further simplified as a Type 1 term can be used when M<K and a Type 5 term when M=K, and reduce the resulting equation to
Figure imgf000071_0003
Calculating residual macrofracture populations
The equations for active and static half-macrofracture populations derived above work well for most of the period, in which the macrofracture network is growing. However, they begin to break down as the fracture network approaches saturation, when the probability of fracture deactivation is high, and the mean lifetime of each macrofracture becomes significantly less that the duration of the timestep (i.e.
FN®0). In this situation, the assumption that fracture deactivation rates are constant for the duration of the timestep is no longer valid.
Therefore, if a simulation of the growth of the fracture network through to saturation is wanted, an alternative method of calculating the half-macrofracture populations in these final timesteps must be found. This can be done by treating the active fracture population as an equilibrium (or“residual”) population, where the rate of fracture nucleation exactly balances the rate of fracture deactivation. Thus, the equilibrium rate of half-macrofracture nucleation can be calculated and used to determine the increase in the static half-macrofracture density during the timestep. The size distribution of the equilibrium half-macrofracture population can also be calculated and used to determine the increase in mean linear half-macrofracture density.
This technique is particularly useful when simulating the growth of an anisotropic fracture network, with a secondary fracture set that develops only after the primary fracture set has already stopped growing. In this situation, the deactivation rate of the secondary fracture set will be high throughout its growth, so the residual growth technique is the only technique that can accurately predict its population distribution function.
Calculating the volumetric density of residual active half-macrofractures MF
The residual active half-macrofracture population is defined as the population where the rate of half-macrofracture deactivation equals the rate of half-macrofracture nucleation, so that the total number of active half-macrofractures remains constant. Therefore
Figure imgf000073_0001
where a rMFP3o(tn) is the residual active half-macrofracture volumetric density, and MFP 3o(tn) is the half-macrofracture nucleation rate, at any time tn during a timestep N. Note that the instantaneous probability of macrofracture deactivation FN must be used, rather than the mean probability of macrofracture deactivation FN, as the mean macrofracture lifetime is short compared with the timestep duration. The half-macrofracture nucleation rate MFP30 can be calculated by differentiating [28] with respect to time, then multiplying the resulting equation by the effective cumulative microfracture activation probability O’, and also by a factor of 2 to convert from macrofractures MF to half-macrofractures JV . This gives:
Figure imgf000073_0002
In deriving [113], the assumption is made that the rate of change of the effective cumulative microfracture activation probability Q’ is negligible compared with the rate of macrofracture nucleation, and hence 0’(tn) can be treated as a constant when differentiating MFP3O. However this does not mean that it can be assumed that the change in the effective cumulative microfracture activation probability over the course of the entire timestep is negligible compared with the change in the rate of macrofracture nucleation during the timestep. In other words, it can be assumed that
Figure imgf000073_0003
but it cannot be assumed that
Figure imgf000074_0001
As the macrofracture network approaches saturation and 0’®0, the variation in 0’(tn) when calculating the rate of macrofracture nucleation at any time tn during timestep N must therefore be taken into account. The rate of change of Q’ can be calculated from the number of propagating macrofracture tips and the rate of tip propagation:
Figure imgf000074_0002
Following a similar derivation to that shown above results in
Figure imgf000074_0003
Inserting CIXMF/CIIIJMF into [116] and integrating gives
Figure imgf000074_0004
[1 18] can now be inserted into [113] to derive a differential equation for the half- macrofracture nucleation rate:
Figure imgf000074_0005
[1 19] is difficult to solve, the problem can be simplified by considering two end- member states: one in which the half-macrofracture nucleation rate is controlled predominantly by changes in the microfracture growth rate, and the other in which the half-macrofracture nucleation rate is controlled predominantly by changes in the effective cumulative microfracture activation probability. It can be determined which of these end-member states to use in any given situation by comparing the residual half-macrofracture volumetric densities calculated for the two end member states, and taking the lowest.
Residual half-macrofracture population controlled by microfracture growth rate
In this scenario, it is assumed that the proportional change in the cumulative microfracture activation probability Q’ during the timestep is negligible compared with the change in the rate at which microfractures grow to radius h/2, i.e.
Figure imgf000075_0001
In this case, [113] can simply be used to calculate the half-macrofracture nucleation rate MFP3O throughout the timestep. Combining this with [112] gives an expression for the total volumetric density of residual active half-macrofractures at any time tn during timestep N:
Figure imgf000075_0002
and the total volumetric density of residual active half-macrofractures at the end of timestep N is given by:
Figure imgf000076_0001
Ideally, O’N should be used rather than O’N-I in [122] However, in practice this is difficult because O’N is dependent on MFP3o(0,tN), so first [122] must be used to calculate MFP3o(0,tN), and then MFP3o(0,tN) must be used to update O’N.
Since macrofracture deactivation is a random process governed by an exponential decay curve, the length distribution of the residual active half-macrofracture population can be calculated by inserting the half-macrofracture length I, propagation rate and instantaneous deactivation probability into the exponential fracture activation probability function [45]. Thus the number of residual active half- macrofractures of length IN or greater is given by
Figure imgf000076_0002
Residual half-macrofracture population controlled by cumulative microfracture activation probability In this scenario it is assumed that the change in the rate at which microfractures grow to radius h/2 is negligible compared with the change in cumulative
microfracture activation probability Q’ during the timestep, i.e.
Figure imgf000076_0003
It can therefore be assumed that the rate at which microfractures grow to radius h/2 is constant throughout the timestep, i.e.
Figure imgf000077_0001
[1 19] can therefore be approximated as
Figure imgf000077_0002
[126] can be solved to give
Figure imgf000077_0003
and therefore the total volumetric density of residual active half-macrofractures at the end of timestep N is given by:
Figure imgf000077_0004
Note that in this case GN can be substituted for GN-I in [128], since G is not dependent on MFP 3O(0,ΪN)·
Calculating the volumetric density of residual static half-macrofractures By definition, the population of active fractures will be small relative to the increase in the static fracture density during the residual growth stage. The active fractures can therefore generally be ignored, and the evolution of the fracture population can be simulated by just calculating the increase in the static half-macrofracture density,
Figure imgf000078_0001
each timestep N.
As previously, to calculate the increase in static half-macrofracture density during a timestep, the active half-macrofracture density is multiplied by the probability of macrofracture deactivation, and the result is integrated throughout the timestep. However this calculation is easier for the residual fracture population, since it is assumed that all the active half-macrofractures have nucleated during the current timestep, so it is not necessary to iterate through previous timesteps of fracture nucleation. As with active residual fractures, there will be two forms of the equations for residual volumetric density of static macrofractures: one form describing residual fracture populations controlled by microfracture growth rate, and one form describing residual fracture populations controlled by the cumulative microfracture activation probability. The form that gives the smallest increase in fracture density will indicate which of these two factors is actually controlling the rate of fracture nucleation at any time.
Residual half-macrofracture population controlled by microfracture growth rate
If the macrofracture nucleation rate is controlled predominantly by changes in the microfracture growth rate, it can be assumed that the cumulative half-macrofracture activation probability Q’ remains constant throughout the timestep. The increase in total static half-macrofracture density during a timestep N can therefore be calculated by inserting [112] and [1 13] into [76]:
Figure imgf000078_0002
Figure imgf000079_0001
This can be multiplied by the residual half-macrofracture length distribution function [123] to derive an expression for the increase in the volumetric density of residual static half-macrofractures of length IN or greater during a timestep N:
Figure imgf000079_0002
Residual half-macrofracture population controlled by cumulative microfracture activation probability
If the macrofracture nucleation rate is controlled predominantly by changes in the cumulative microfracture activation probability, the integral in [129] is more complex as the change in the cumulative half-macrofracture activation probability Q’ throughout the timestep must be taken into account. However, [127] can be inserted into the equation
Figure imgf000079_0003
and the result can be integrated to obtain an expression for the total incremental increase in the volumetric density of residual static half-macrofractures
AS I-MFP 3o(0, N). This gives:
[132] GMI¾O (0’ N)
Figure imgf000080_0001
Again, this can be multiplied by the residual half-macrofracture length distribution function [123] to derive an expression for the incremental increase in the volumetric density of residual static half-macrofractures of length IN or greater during timestep
N:
Figure imgf000080_0002
Calculating the mean linear density of residual active half-macrofractures
The mean linear density of residual half-macrofractures can be calculated in the same way as above: by differentiating the appropriate MFP30 functions with respect to I to determine the number of half-macrofractures with length between I and l+dl, multiplying this by Ih to obtain their total area, and then integrating the resulting expression from IN to the maximum possible half-macrofracture half-length ax. However this is much simpler for residual half-macrofractures, because the cumulative density distribution functions for residual half-macrofractures all follow the form of [123]:
[134] rMFPso O^JVi tjy) = rMF¾Q (0 tjv)
Figure imgf000081_0001
The differentiation and integration on the length-dependent exponential term can therefore be carried out independently of the expression for total residual half- macrofracture volumetric density:
Figure imgf000081_0002
Thus the total mean linear density of all residual active half-macrofractures at the end of timestep N is given by
Figure imgf000081_0003
and the mean linear density of residual active half-macrofractures of length IN or greater at the end of timestep N is given by
Figure imgf000081_0004
where a rMFP3o(0,tN) is given by [122] or [128] as appropriate.
Calculating the mean linear density of residual static half-macrofractures The same method can be used to calculate the incremental increase in the mean linear density of residual static half-macrofractures during a timestep N, AS,-MFP (N). Thus the incremental increase in the total mean linear density of all residual static half-macrofractures during timestep N is given by
Figure imgf000082_0001
and the incremental increase in the mean linear density of residual static half- macrofractures of length IN or greater during timestep N is given by
Figure imgf000082_0002
where D5 GMER3O(0,N) is given by [129] or [132] as appropriate.
Elastic moduli and stress
One of the key features of this method is that it takes into account the effect that the growing fracture populations will have on the internal elastic stress field within the layer (the in situ stress). This is important because the in situ stress drives fracture propagation, and controls the rate of fracture propagation through the fracture driving stress Od. However the growth of the fractures will also affect the in situ stress because fracture displacement accommodates some of the horizontal strain that would otherwise be accommodated elastically in the rock mass, either locally (in the form of a stress shadow) or distributed throughout the layer. This creates a feedback mechanism that is likely to be an important control on the final fracture density and distribution.
The boundary conditions for tectonic deformation can typically be defined in terms of a constant effective vertical stress ov’, and a bimodal horizontal strain applied at specified rates έ Ghίh and έpGhqc, in a specified orientation and for a specified duration. In a perfectly elastic rock mass, this information, combined with the bulk rock elastic moduli, is sufficient to calculate the driving stress on each fracture set, at any time t. If the horizontal strain is applied at a constant rate, the fracture driving stress will increase linearly with time. In practice, however, a rock mass subject to strain over geological timescales will often experience time-dependent (viscoelastic) deformation processes such as creep, grain sliding and pressure solution, which gradually convert elastic strain into inelastic strain. In the initial stages of deformation, these“strain relaxation” processes will be slow and, thus, the elastic strain in the rock mass will increase in proportion to the applied strain rate. However, as the elastic strain increases, the strain relaxation processes will accelerate until an equilibrium is reached between the applied strain rate and the rate of strain relaxation. Thus, the deformation will be characterised by an early stage with increasing fracture driving stress, followed by a late stage with constant fracture driving stress. The ultimate fracture geometry will depend on the rate of fracture growth relative to the rate of strain relaxation:
fractures networks may grow rapidly in the early stage, with later strain
accommodated by increasing displacement on the existing fractures, or they may grow slowly in the late stage. For timesteps in which the fracture driving stress is increasing, it is useful to define a “weighted mean” driving stress, such that
Figure imgf000083_0001
where OdN-i is the fracture driving stress at the end of the previous timestep N-1 , and dod/dt is the rate of change of fracture driving stress during the current timestep N. Since the rate of fracture propagation is proportional to the fracture driving stress to the power of the subcritical fracture propagation index b, the weighted mean driving stress odN can be substituted for OdN, to calculate the overall growth of the fracture populations throughout the timestep. The bulk rock elastic moduli take account of the elastic strain accommodated on the fractures as well as the elastic strain in the rock mass. Therefore, they may themselves change through time as the fracture network grows. The exact form of the bulk rock elastic moduli will depend on the distribution of stress throughout the layer: whether there are fracture stress shadows or the stress is evenly distributed.
Elastic moduli in the stress shadow scenario
In this scenario, elastic displacement on the fractures accommodates all the elastic strain from the rock mass in a zone immediately surrounding the fracture. This zone (the stress shadow) will therefore be characterised by a local in situ stress anomaly, while the stress and strain in the rock mass further away from the fractures will be uniform, reflecting the regional applied strain and the elastic properties of the rock mass. Since only fractures that lie outside the stress shadows can propagate, the horizontal effective stresses driving the propagation of these fractures will be independent of the fracture population, and can be calculated from the effective vertical stress ov’, the horizontal elastic strains Shmin and Shmax, and the elastic properties of the intact rock Er (Young’s Modulus) and vr (Poisson’s ratio):
Figure imgf000084_0001
where ehmm(t) and Shmax(t) represent the instantaneous elastic horizontal strains at time t, which will be a function of the applied strain rates and the strain relaxation rate.
Note that while a negative (extensional) strain is required to generate a driving stress for Mode 1 fractures 101 , it is possible to have a driving stress for Mode 2 fractures 102 even with no applied external horizontal strain (i.e. when
ehmin=ehmax=0), if there is a sufficient differential stress (ov’-Ohmm’). This suggests that Mode 2 shear fractures 102 can develop purely in response to the lithostatic load, and may correspond to the inferred compaction-related deformation bands observed in some sandstones, chalks and other lithologies. It is likely, however, that these early fractures will heal over geological time, and will not be reactivated during tectonic deformation. Equation [141] can account for the resulting compaction- related strain through the initial compaction factor u. u can vary between u=0, which represents an elastic equilibrium, and u=1 , which represents a viscoelastic equilibrium where ohmin’ = ahmax’ = sn’ when ehmm= ehmax = 0.
Elastic moduli in the evenly distributed stress scenario
In this scenario, elastic displacement on the fractures accommodates some of the applied external horizontal strain, and the remainder of the strain is evenly distributed throughout the rock mass. No local stress anomalies develop around the fractures.
In this case, the fractures will introduce an anisotropy in the elastic response of the bulk rock to an applied external strain. If it is assumed that the fracture sets are orthogonal to the principal horizontal stress directions, the orthotropic form of Hooke’s Law can be used to relate bulk rock elastic strain to in situ stress:
Figure imgf000085_0001
Based on the fracture strain expressions [38], the expressions for the bulk rock elastic constants in [142] can be derived:
[143]
Figure imgf000085_0002
Figure imgf000085_0003
and
Figure imgf000086_0001
where hmin MFP32 represents the mean linear density of all fractures (active and static) striking perpendicular to the Chmin orientation, hma)WP32 represents the mean linear density of all fractures striking perpendicular to the Shmax orientation, and Mhh, MVh, Mhv and Mw are a set of Fracture Mode factors defined as follows:
• For vertical Mode 1 dilatant fractures 101 striking perpendicular to hmin
[145]
Figure imgf000086_0002
• For vertical Mode 1 dilatant fractures 101 striking perpendicular to hmax
[146]
Figure imgf000086_0003
• For normal Mode 2 shear fractures 102 with dip w and friction coefficient pfr
[147] Mhh = cos w (sin w
Mvh =— cos w (sin
Mhv =— cos w (sin
Figure imgf000086_0004
Figure imgf000087_0001
The constitutive relationships can be simplified considerably by defining fracture multipliers for the Young’s Modulus and Poisson’s ratio, E and v respectively:
Figure imgf000087_0002
Now, the bulk rock constitutive relationship can be written as
Figure imgf000087_0003
Note that the bulk moduli defined in [143] and [144], and the fracture multipliers defined in [148] and [149], do not take account of the microfractures, which are assumed to accommodate a negligible amount of strain in comparison to the layer- bound macrofractures. However, in situations where the microfractures are significant, it is possible to expand [143], [144], [148] and [149] to include additional terms representing elastic strain accommodated on the microfractures. These terms will be based on the total fracture porosity for Mode 1 microfractures or the total volumetric heave for Mode 2 microfractures, which must be calculated numerically using a similar procedure to that used to calculate the mean linear microfracture density MFP32.
Application of the method to geological formations
In order to apply this method to fractured layers in real geological formations, a workflow must be described that includes all the elements of the above-described calculations or similar calculations in an appropriate order. This will allow for the prediction of the fracture populations in a specified layer as a function of the mechanical properties and deformation history (effective vertical stress, horizontal strain rate, and duration of deformation).
As mentioned above, the method uses a combined implicit and explicit fracture model. The implicit fracture model is used for simulating the processes of fracture nucleation, propagation and interaction, based on geomechanical principles and representing fracture populations statistically by functions or values. The explicit fracture model, on the other hand, is used for simulating fracture propagation, interactions and intersections between different fractures in a Discrete Fracture Network (DFN), in which a fracture is represented as an individual object, based on data from the implicit fracture model, such as fracture nucleation, propagation and growth rates.
The calculations described so far assume that the layer is planar and flat, and that the mechanical properties and deformational history are homogeneous and extend laterally for an infinite distance. By contrast, real geological layers tend to have complex geometries and laterally variable mechanical properties and deformation histories - consider for example the variation in total strain between a fold hinge and a fold limb, or the complex fracture patterns that can develop in a fault damage zone. For this reason, static geological models typically discretise the 3D geometry into a grid of discrete (usually hexahedral) gridblocks. The mechanical properties and in situ stress are assumed to be uniform within each gridblock, but can vary between gridblocks.
This complexity and variability can be incorporated into the present fracture models by running the implicit components of the method independently for each gridblock, using the mechanical properties and deformation history defined for that gridblock. This will allow for an easy comparison of the fracture density and geometry predicted in different parts of the geological structure - e.g. the predicted mean fracture spacing and length in a fold hinge and a fold limb, or in the footwall and hanging wall of a fault can be compared. However, to obtain accurate results, there are some constraints on the gridblock geometry: · The gridblocks should extend vertically across the entire thickness of the fractured layer - i.e. the thickness of the layer should correspond to the height of one gridblock.
The lengths and widths of the gridblocks should be greater than the typical maximum length of the macrofractures, so that the majority of macrofractures that nucleate within a gridblock will remain entirely within that gridblock.
The thickness and depth of the layer should not vary by too much within each gridblock - i.e. they should be relatively flat-lying (or at least should have been at the time of deformation), and should not have significant taper.
It may be necessary to upscale the gridblocks in the original static geomodel to meet these criteria. The implicit fracture data for each gridblock can be combined to generate a single, integrated explicit DFN spanning the entire geological model. This allows for the generation of large, geomechanically consistent DFNs that reflect the variability in geometry, mechanical properties and in situ stress across a large geological structure. Fig. 8 shows the overall process steps of the implicit fracture modelling procedure according to an embodiment of the invention. This procedure is used for generating implicit fracture population data for each gridblock in the static geomodel, in the form of cumulative fracture density functions MFP3o(r), MFP32(r), MFP3O(I) and MFP32(I). The procedure can also be used to generate other useful data described in the previous sections, such as the fracture aperture and porosity, total fracture strain, mean fracture length, and fracture connectivity. This data can be broken down by fracture set if required.
Step 801 : The implicit fracture modelling procedure is initiated by retrieving data describing the geomechanical model. Such data includes, for instance, vertical and lateral dimensions of each of the geological layers, mechanical property data for the geological layers and stress and strain history data.
Step 802: The model is discretised into a number of gridblocks.
Steps 803 and 808: The remaining steps 804-807 are worked through for all gridblocks, one by one.
Step 804: The initial stress state is calculated, based on the depth of burial, on the fluid pressure and on the degree of initial compaction. Equations [141]-[151] can be used for this calculation.
Step 805: The evolution of the total active and static fracture densities MFP3O(0), MFP 32(0), MFP3O(0) and MFP32(0) for each fracture set are calculated, looping through the timesteps.
Step 806: The cumulative macrofracture populations are calculated, looping through the timesteps.
Step 807: The cumulative microfracture populations are calculated, looping through the timesteps. Steps 805-807 are described in more detail in the following in relation to Figs. 9-1 1 , respectively.
Fig. 9 shows in more detail the process steps to calculate total fracture densities according to an embodiment of the invention. Thus, Fig. 9 is a breakdown of step 805.
Steps 901 and 910: The individual steps 902-909 are worked through for all timesteps or until the fracture network is complete, one timestep at a time.
Step 902: The stress state is recalculated at the start of each timestep. This recalculation is based on the vertical stress, on the fluid pressure, on the total horizontal strain and on the strain relaxation. It also takes account of the impact of fracture growth on the mechanical properties, using fracture populations from the previous timestep. Like in step 804, in which the initial stress state is calculated, Equations [141]-[151] can be used for this calculation.
Step 903: The optimal timestep duration is estimated for each timestep. Thus, the duration of each timestep is determined dynamically to ensure that the maximum increase in macrofracture stress shadow volume remains below a specified limit; this ensures that the used approximation that fracture deactivation rates are constant throughout the timestep is valid. The estimation is can be done by calculating, for each fracture set, the time taken for the macrofracture density to grow by a specified amount (typically between 0.2% and 1 %) if there is no fracture deactivation and then choosing the smallest of the calculated time values.
Step 904: Based on the recalculated stress state, on the estimated optimal timestep duration, on the applied strain rate and on the fracture set orientation, the mean fracture driving stress, i.e. the local stress acting on each fracture set (Equation [140]), the growth factor parameters y and G (Equations [7]-[12]) and the maximum macrofracture propagation distance OMFA (Equations [19]-[29]) are calculated.
Step 905: The rate of macrofracture deactivation is calculated. Macrofracture deactivation can occur by two mechanisms: Interaction with the stress shadow of other parallel macrofractures or intersection with other non-parallel macrofractures.
In both cases, the probability that a specific macrofracture will be deactivated during a given timestep is a function of the fracture growth rate and the fracture populations (at the end of the previous timestep). Equations [50]-[65] can be used for the calculation of the macrofracture deactivation probabilities F and F.
Step 906: The increase in the total active and static macrofracture densities MFP3O(0) and MFP 32(0) are calculated analytically for each fracture set. These calculations are based on the timestep duration, on the driving stress, on the growth factor, and on the macrofracture deactivation probability for each fracture set. Equations [72]-[11 1] can be used for these calculations. If the macrofracture deactivation probability is high, it may be necessary to use the formulae for residual macrofracture populations (Equations [1 12]-[139]).
Step 907: The rate of microfracture deactivation is calculated. Microfracture deactivation can occur when a microfracture falls into the stress shadow of a macrofracture. Thus, it is a function of the growth of the macrofracture population in the timestep and it can only be calculated after the calculation of the macrofracture growth. Equations [45]-[49] can be used for the calculation of the microfracture deactivation probabilities q and Q.
Step 908: The increase in the total active and static microfracture densities mrR3o(0) and mrR 32(0) are calculated analytically and numerically, respectively, for each fracture set. These calculations are based on the timestep duration, on the driving stress, on the growth factor, and on the microfracture deactivation probability for each fracture set. Equations [66]-[71] can be used for these calculations.
Step 909: The calculation continues until the fracture network is complete, i.e. until all fracture sets have stopped growing or the specified maximum time or timestep limit has been reached. A fracture set is assumed to have stopped growing, if the density MFP3O(0) is below a specified value, or if the total exclusion zone volume approaches 100%. Fig. 10 shows in more detail the process steps to calculate cumulative
macrofracture density distribution functions according to an embodiment of the invention. Thus, Fig. 10 is a breakdown of step 806.
Steps 1001 and 1004: The individual steps 1002 and 1003 are worked through for all timesteps, one by one.
Step 1002: Dynamic fracture data are retrieved for each timestep. These data include the timestep duration, the fracture driving stress, the growth factor and the macrofracture deactivation probability for each fracture set, all of which were calculated and saved in step 805.
Step 1003: The increase in the cumulative active and passive macrofracture densities as a function of length MFP3O(I) and MFP32(I) are calculated analytically for a range of specified half-lengths I. These calculations are based on the timestep duration, on the driving stress, on the growth factor, and on the macrofracture deactivation probability for each fracture set. Equations [72]-[11 1] can be used for these calculations. If the macrofracture deactivation probability is high, it may be necessary to use the formulae for residual macrofracture populations (Equations [1 12]-[139]). The combination of results for all specified half-macrofracture lengths I define the cumulative macrofracture length distribution functions.
Fig. 1 1 shows in more detail the process steps to calculate cumulative microfracture densities according to an embodiment of the invention. Thus, Fig. 1 1 is a breakdown of step 807.
Steps 1101 and 1104: The individual steps 1102 and 1 103 are worked through for all timesteps, one by one.
Step 1 102: Dynamic fracture data are retrieved for each timestep. These data include the timestep duration, the fracture driving stress, the growth factor and the microfracture deactivation probability for each fracture set, all of which were calculated and saved in step 805. Step 1 103: The increase in the cumulative active and passive microfracture densities as a function of radius uFP3o(r) and ^32( are calculated for a range of specified fracture radii r. These calculations are based on the timestep duration, on the driving stress, on the growth factor, and on the microfracture deactivation probability for each fracture set. ^30( is calculated analytically and ^32( is calculated by numerical integration between the specified radii r. Equations [66]-[71] can be used for these calculations. The combination of results for all specified radii r define the cumulative microfracture size distribution functions.
Fig. 12 shows the overall process steps of the explicit fracture modelling procedure according to an embodiment of the invention. This procedure is used for generating the necessary data for construction of the explicit DFN model, based on the calculated implicit fracture data calculated for each of the gridblocks in the geomodel as described above. The DFN model can be populated either with both
microfractures and macrofractures or with macrofractures only. In the latter case, seed macrofractures are added to the DFN according to the macrofracture nucleation rates calculated with the implicit fracture data. This will reduce the total number of fractures in the DFN while retaining all the long fractures that are likely to act as major fluid flow pathways, thus making it more suitable as input for fluid flow modelling.
The procedure uses the timesteps defined for each gridblock when calculating the implicit fracture data. It loops through all of the timesteps in all of the gridblocks in the geomodel in chronological order. Thus, it can use the driving stress, nucleation rate and propagation rate data calculated for each fracture set in each timestep in the implicit calculation phase to control the nucleation and growth of the fractures in the explicit DFN.
Step 1201 : The explicit fracture modelling procedure is initiated by retrieving implicit fracture data and relevant dynamic fracture data for each fracture set and each timestep from the calculated implicit fracture model. The dynamic fracture data includes, for instance, mean fracture growing stress, growth factor and maximum fracture propagation distance. Step 1202: A sorted list of calculation elements in the form of“gridblock-timesteps” is created. Each such gridblock-timestep consists of a specific timestep of a specific gridblock. Since the timestep durations are calculated individually for each gridblock, they will vary between gridblocks. Therefore, a list of all such gridblock-timesteps is created and sorted in chronological order by the end times of each gridblock- timestep.
Steps 1203 and 1211 : The remaining steps 1204-1210 are worked through for all gridblock-timesteps, one by one.
Steps 1204 and 1210: The remaining steps 1205-1209 are worked through for all fracture sets, one by one.
Step 1205: If microfractures are to be included in the explicit DFN, steps 1206 and 1207 are performed and step 1208 is skipped. If not, steps 1206 and 1207 are skipped and step 1208 is performed.
Step 1206: New microfractures are added.
Step 1207: The growth of each microfracture (newly added ones as well as already existing ones) is calculated, including a possible transformation into a macrofracture.
Step 1208: New macrofractures are added directly into the model.
Step 1209: The growth of each macrofracture is calculated, taking into account possible interactions with parallel macrofractures and intersections with non-parallel macrofractures as well as possible crossings of gridblock boundaries.
Steps 1206-1209 are described in more detail in the following in relation to Figs. 13- 16, respectively.
Fig. 13 shows in more detail the process steps to add new microfractures according to an embodiment of the invention. Thus, Fig. 13 is a breakdown of step 1206. Step 1301 : It is calculated how frequently microfractures will grow to a specified minimum radius. Since the total microfracture density is infinite for a power law distribution, a minimum radius limit for the microfractures in the DFN must be specified. Equation [14] can be used to determine how many microfractures will reach the minimum radius during a given timestep. This calculation is based on the microfracture growth rate, on the cumulative microfracture density distribution functions and on the gridblock volume.
Step 1302: The current time is set to the start of the timestep being processed. If the stress is not constant, a weighted time can be used.
Step 1303: A random location for a new microfracture is picked within the gridblock.
Step 1304: It is checked whether the picked location lies within a stress shadow of a macrofracture of the same fracture set.
Step 1305: If the picked location lies within such a stress shadow, no new microfracture will be added to the model.
Step 1306: If the picked location does not lie within such a stress shadow, a new microfracture is created at the picked location. The nucleation time of the new microfracture is set to be the current time (or alternatively the weighted time) and the radius of the new microfracture is set to be the specified minimum radius.
Step 1307: The current time is incremented according to the calculated frequency at which microfractures will grow to a specified minimum radius. Alternatively, the weighted time is used.
Step 1308: Steps 1303-1307 are repeated until the current time (or alternatively the weighted time) has reached the end of the timestep being processed.
Fig. 14 shows in more detail the process steps to calculate microfracture growth according to an embodiment of the invention. Thus, Fig. 14 is a breakdown of step 1207. Steps 1401 and 1410: The individual steps 1402-1409 are worked through for all microfractures, newly added as well as already existing, one by one. Step 1402: It is checked whether the microfracture lies within a stress shadow of a macrofracture of the same fracture set. The microfracture can have fallen into such a stress shadow since the previous timestep.
Step 1403: If the microfracture lies within such a stress shadow, the microfracture is deactivated.
Step 1404: If the microfracture is not active, steps 1405-1409 are skipped
Step 1405: The microfracture radius is incremented. The new microfracture radius is a function of the current radius, of the microfracture growth rate and of the growth time, the growth rate being determined by equation [8] and the fracture driving stress for the timestep being processed. For microfractures already present at the start of the timestep, the growth time is the timestep duration, whereas, for microfractures added during the timestep (in step 1306), the growth time is the time passed from the nucleation time to the end of the timestep.
Step 1406: It is checked whether the microfracture has grown into a macrofracture, i.e. if the new radius is greater than half the layer thickness. Step 1407: If the microfracture has not grown into a macrofracture, steps 1408 and 1409 are skipped.
Step 1408: A new macrofracture of zero length is generated at the location of the microfracture.
Step 1409: The microfracture is deactivated. Fig. 15 shows in more detail the process steps to add new macrofractures directly according to an embodiment of the invention. Thus, Fig. 15 is a breakdown of step 1208.
Step 1501 : It is calculated how frequently macrofractures will nucleate from growing microfractures, based on data from the implicit fracture model. New macrofractures will nucleate whenever a microfracture grows to a radius of half the layer thickness, cf. step 1406 above. Thus, the macrofracture nucleation rate is a function of the microfracture growth rate, of the cumulative microfracture density distribution functions, of the gridblock volume, and of the layer thickness. If microfractures are not included in the DFN, new seed macrofractures can be added directly, at a rate determined by equation [28], with I_N=0.
Step 1502: The current time is set to the start of the timestep being processed. If the stress is not constant, a weighted time can be used.
Step 1503: A random location for a new macrofracture is picked within the gridblock.
Step 1504: It is checked whether the picked location lies within an exclusion zone of a macrofracture of the same fracture set.
Step 1505: If the picked location lies within such an exclusion zone, no new macrofracture will be added to the model.
Step 1506: If the picked location does not lie within such an exclusion zone, a new macrofracture is created at the picked location. The nucleation time of the new macrofracture is set to be the current time (or alternatively the weighted time) and the length of the new macrofracture is set to be zero.
Step 1507: The current time is incremented according to the calculated frequency at which macrofractures will nucleate. Alternatively, the weighted time is used.
Step 1508: Steps 1503-1507 are repeated until the current time (or alternatively the weighted time) has reached the end of the timestep being processed. Figs. 16a and 16b show in more detail the process steps to calculate macrofracture growth according to an embodiment of the invention. Thus, Figs. 16a and 16b are a breakdown of step 1209.
Steps 1601 and 1619: The individual steps 1602-1618 are worked through for all macrofracture tips, of newly added macrofractures as well as of already existing macrofractures, one by one.
Step 1602: If the macrofracture tip is not active, steps 1603-1618 are skipped.
Step 1603: The maximum propagation distance of the macrofracture tip is calculated. All macrofracture tips within a specific fracture set propagate at the same rate. The propagation rate can be calculated using equation [20], substituting e for L/2. Thus, the maximum propagation distance for a given macrofracture tip within a given timestep is the propagation rate times the propagation time within the timestep. For macrofractures already present at the start of the timestep, the propagation time is the timestep duration, whereas, for macrofractures added during the timestep (in step 1506), the propagation time is the time passed from the nucleation time to the end of the timestep.
Step 1604: If the propagating macrofracture tip will interact with a stress shadow of another macrofracture of the same fracture set if growing by the maximum propagation distance, steps 1605-1607 are performed. If not, these three steps are skipped. Such interaction occurs if the stress shadows of the two macrofractures overlap.
Step 1605: The maximum propagation distance is reduced to the distance, at which stress shadow interaction occurs.
Step 1606: Both of the interacting macrofracture tips are deactivated.
Step 1607: If required, a relay segment is added to link the two interacting macrofractures. Step 1608: If the propagating macrofracture tip will intersect a macrofracture of another fracture set if growing by the maximum propagation distance, steps 1609 and 1610 are performed. If not, these two steps are skipped.
Step 1609: The maximum propagation distance is reduced to the distance, at which intersection occurs.
Step 1610: The propagating macrofracture tip is deactivated.
Step 161 1 : If the propagating macrofracture tip will cross a gridblock boundary if growing by the maximum propagation distance, steps 1612-1617 are performed and step 1618 is skipped. If not, steps 1612-1617 are skipped and step 1618 is performed.
Step 1612: The real time, at which the propagating macrofracture tip will cross the gridblock boundary, is calculated.
Step 1613: The maximum propagation distance is reduced to the distance, at which the gridblock boundary is reached.
Step 1614: The propagating macrofracture tip is extended by the maximum propagation distance.
Step 1615: The propagating macrofracture tip is deactivated.
Step 1616: A new macrofracture of zero length is added in the adjacent gridblock at the entering position. The nucleation time is set to the time calculated in step 1612. One of the macrofracture tips of the new macrofracture is connected to the deactivated macrofracture tip at the gridblock boundary. The status of the other macrofracture tip of the new macrofracture is set to active.
Step 1617: If fracture propagation has already been calculated up to this time in the adjacent gridblock, the active macrofracture tip of the new macrofracture is extended in the appropriate direction of propagation into the adjacent gridblock until the position corresponding to the latest calculation time of the adjacent gridblock, checking for stress shadow interaction and intersection as described above in steps 1604-1610. It should be noted that the orientation and propagation rate of macrofractures as well as the timestep durations may be different in the adjacent gridblock. Note that if the fracture orientations in the two adjacent gridblocks are convergent, it will be necessary to extend the fracture along the gridblock boundary to avoid crossing backwards and forwards across the gridblock boundary in an infinite loop. If fracture propagation has not yet been calculated up to this time in the adjacent gridblock, the new macrofracture is left at zero length to be calculated in the next timestep for the adjacent gridblock.
Step 1618: The propagating macrofracture tip is extended by the maximum propagation distance.
Fig. 17 shows the process steps to construct an explicit DFN according to an embodiment of the invention. The macrofractures generated by the above-described modelling procedures are segmented. This allows the DFN to replicate complex macrofracture geometries, for example macrofractures with variable height or curved macrofractures that follow local variations in the in situ stress orientation (as is often observed in outcrop). Once the simulation of macrofracture propagation in all the gridblocks has been finished, the connected planar segments must therefore be recombined into single, generally non-planar entities representing entire macrofractures.
Step 1701 : All macrofractures, which are connected across gridblock boundaries are combined to create continuous macrofractures that span across multiple gridblocks.
Step 1702: The positions of the corner points of all macrofractures are calculated to ensure accurate bevelling across gridblock boundaries and macrofracture intersections. Such bevelling is needed for ensuring a smooth alignment of the connected macrofracture segments. Fig. 18 is a schematic illustration of bevelling of corner points between two non- parallel macrofracture segments according to an embodiment of the invention. In Fig. 18, P1 and P2 are the two top corners of a first macrofracture segment P propagating towards a gridblock boundary GB in a direction from P1 to P2. Q1 and Q2 are the two top corners of a second macrofracture segment Q propagating away from the gridblock boundary in a direction from Q1 to Q2.
It is clear that the top corners P2 and Q1 of the macrofracture segments P, Q, respectively, are not coincident. Therefore, an intersection point T between the top edges of the macrofracture segments P, Q, respectively, is calculated and used as the corner point for both macrofracture segments P, Q.
If the four top corner points P1 , P2, Q1 and Q2 of the macrofracture segments P, Q are given by (P1x, P1y, P1z), (P2x, P2y, P2z), (Q1x, Q1y, Q1z) and (Q2x, Q2y, Q2z), respectively, in global xyz coordinates, then the following can be defined
[152] APx = P2x— Plx
APy = Ply— Ply
APz = P2z— Plz
AQx = Q2x— Qlx
AQy = Q2y - Qly
AQz = QZz— Qlz and the intersection point T between the top edges of the macrofracture segments P, Q will be given by
Figure imgf000102_0001
A similar formula can be used to calculate the intersection point between the bottom edges of the two macrofracture segments P, Q. LIST OF CONSTANTS USED IN THE CALCULATIONS
The following fundamental geometric and material properties are assumed to be constant in the calculations. Where appropriate, typical values are given for consolidated sedimentary rocks at c.2km burial depth.
Figure imgf000103_0001
LIST OF SYMBOLS USED IN THE CALCULATIONS
The following symbols represent derived variables that are parameters of or output from the calculations described in this document. In some cases, modifiers
(subscripts or superscripts) are applied to constrain the scope of the variables (for example, so they apply only to a specific fracture type or set, or a specific timestep). Where these modifiers are absent, the variables can be assumed to apply generally or in situations in which the scope is not defined; e.g. if one fracture type or set is considered, if the timestep is undefined.
Figure imgf000104_0001
Figure imgf000105_0001
Figure imgf000106_0001
LIST OF REFERENCES USED IN THE DRAWINGS
101 Mode 1 dilatant fracture
102 Mode 2 shear fracture
201 Stress shadow around first macrofracture
301 Stress shadow around the second macrofracture
302 Exclusion zone around second macrofracture
303 Stress shadow interaction box
A Macrofracture stress shadow interaction
B Macrofracture intersection
J Timestep of nucleation of half-macrofracture with specific length
GB Gridblock boundary
K Timestep of fracture deactivation
I Length of half-macrofracture
IN Specific length of half-macrofracture at time tn
LSSIB Length of stress shadow interaction box
LT Lateral macrofracture tip
LTA Lateral macrofracture tip of first macrofracture
LTB Lateral macrofracture tip of second macrofracture
LTC Lateral macrofracture tip of third macrofracture
M Timestep of fracture nucleation
MF Macrofracture
MFA First macrofracture
MFB Second macrofracture
MFc Third macrofracture
P First macrofracture segment
P1 First top corner of first macrofracture segment
P2 Second top corner of first macrofracture segment
Q Second macrofracture segment
Q1 First top corner of second macrofracture segment
Q2 Second top corner of second macrofracture segment
t Time
T Intersection point between first and second macrofracture segments W Mean stress shadow width mR Microfracture
dhmax Maximum principal horizontal strain CJhmin Minimum principal horizontal strain sn’ Vertical stress
w Dip angle

Claims

1. A method for modelling and simulating a fractured geological structure using a dynamic geomechanical model, which is discretised into geological layers, timesteps and gridblocks and comprises an implicit fracture model and an explicit discrete fracture network model, wherein the method comprises an implicit fracture modelling procedure for simulating the processes of fracture nucleation, propagation and interaction, based on fundamental geomechanical principles, to generate an implicit fracture model, in which model a fracture population is represented statistically by functions or values, and an explicit fracture modelling procedure for simulating the process of fracture propagation, interactions between parallel fractures and
intersections of non-parallel fractures in an explicit discrete fracture network model, in which a fracture is represented as an individual object, based on data from the implicit fracture model relating to fracture nucleation and propagation (such as rates of fracture nucleation and growth), as well as incorporating the explicit discrete fracture network model into a geological simulator.
2. The method according to claim 1 , wherein the dynamic geomechanical model comprises one or more microfractures contained entirely within a geological layer.
3. The method according to claim 2, wherein the one or more microfractures are circular and the characterising dimension of each of the microfractures is a radius.
4. The method according to any of the preceding claims, wherein the dynamic geomechanical model comprises one or more macrofractures, whose upper and lower edges lie on the upper and lower surfaces of a geological layer.
5. The method according to claim 4, wherein the macrofractures are
rectangular or polygonal in shape and the characterising dimension of each of the macrofractures is a length in a direction parallel to the layer by the surfaces of which the macrofracture is limited.
6. The method according to any of the preceding claims, wherein the implicit fracture modelling procedure comprises the steps of a. retrieving data describing a geomechanical model associated with a geologic environment; b. discretising the geomechanical model into a number of gridblocks; c. performing the following steps for each of the gridblocks: i. calculating an initial stress state; repeating each of the following until a specified time or a specified number of timesteps is reached: ii. calculating total fracture densities;
iii. calculating cumulative macrofracture density distribution functions for each fracture set;
IV. calculating cumulative microfracture density distribution functions for each fracture set.
7. The method according to claim 6, wherein the data describing a
geomechanical model comprises the following: a. vertical dimension and lateral dimensions of each of the geological layers;
b. mechanical property data for the geological layers; c. stress and strain history data.
8. The method according to claim 6 or 7, wherein the step of calculating total fracture densities is performed by repeating the following steps for each timestep until the fracture network of the gridblock is complete: a. recalculating the stress state;
b. estimating an optimal timestep duration;
c. calculating and saving a mean fracture driving stress, a growth factor and a maximum macrofracture propagation distance for each fracture set;
d. calculating a rate of macrofracture deactivation for each fracture set; e. calculating the total increase in active and static macrofracture
density distribution functions for each fracture set;
f. calculating a rate of microfracture deactivation for each fracture set; g. calculating the total increase in active and static microfracture density distribution functions for each fracture set;
h. checking if the fracture network of the gridblock is complete, i.e. if all macrofractures in all fracture sets have stopped growing or if the specified time or the specified number of timesteps has been reached.
9. The method according to claim 8, wherein the estimation of the optimal timestep duration is done by calculating, for each fracture set, the time taken for the macrofracture density to grow by a specified amount (typically between 0.2% and 1%) if there is no fracture deactivation and then choosing the smallest of the calculated time values.
10. The method according to claim 8 or 9, wherein the calculation of the mean fracture driving stress, the growth factor and the maximum macrofracture propagation distance is based on the estimated optimal timestep duration.
1 1. The method according to any of claims 8-10, wherein the calculation of the rate of macrofracture deactivation is based on the calculated growth factor, the calculated maximum macrofracture propagation distance and the total macrofracture populations as calculated at the end of the previous timestep.
12. The method according to any of claims 8-1 1 , wherein the calculation of the total increase in active and static macrofracture populations during the timestep is based on the estimated optimal timestep duration, the calculated mean fracture driving stress, the calculated rate of fracture deactivation, the calculated growth factor and the calculated maximum macrofracture propagation distance for the timestep.
13. The method according to any of claims 6-12, wherein the step of calculating the cumulative macrofracture density distribution functions for each fracture set is performed by repeating the following steps for each timestep: a. retrieving relevant dynamic fracture data (such as mean fracture driving stress, growth factor, maximum macrofracture propagation distance and macrofracture deactivation rate) for the relevant timestep;
b. calculating the increase in the cumulative macrofracture density distribution functions for a range of specified half-macrofracture lengths.
14. The method according to any of claim 6-13, wherein the step of calculating the cumulative microfracture density distribution functions for each fracture set is performed by repeating the following steps for each timestep: a. retrieving relevant dynamic fracture data (such as mean fracture driving stress, growth factor and microfracture deactivation rate) for the relevant timestep;
b. calculating the increase in the cumulative microfracture density
distribution functions for a range of specified fracture radii.
15. The method according to any of the preceding claims, wherein the explicit fracture modelling procedure comprises the steps of a. retrieving implicit fracture data and relevant dynamic fracture data (such as mean fracture driving stress, growth factor and maximum macrofracture propagation distance) for each fracture set and each timestep from the calculated implicit fracture model; b. creating a sorted list of gridblock-timesteps, each consisting of a specific timestep of a specific gridblock, by generating a list of all timesteps in all gridblock and sorting it by the end times of the timesteps; c. performing the following steps for each fracture set in each gridblock- timestep: i. if microfractures are to be included in the explicit discrete fracture network:
1. adding new microfractures;
2. calculating microfracture growth; ii. if microfractures are not to be included:
1. adding new macrofractures directly iii. calculating macrofracture growth; d. constructing the explicit discrete fracture network for the whole grid.
16. The method according to claim 15, wherein the step of adding new
microfractures is performed by: a. calculating how frequently microfractures will grow to a specified minimum radius, based on implicit fracture data (such as cumulative density distribution functions) and relevant dynamic fracture data (such as mean fracture driving stress and growth factor) for the relevant fracture set and timestep; b. setting current time to start of the timestep; c. repeating the following steps until the end of the timestep has been reached: i. picking a random location for a new microfracture;
ii. checking if the location for the new microfracture lies within a stress shadow of a macrofracture of the same fracture set; iii. if the location does not lie within such a stress shadow:
1. adding a new microfracture to the explicit discrete fracture network at the picked location, setting the status of the new microfracture to active;
iv. incrementing the current time by an amount of time
corresponding to the calculated frequency of microfractures growing to the specified minimum radius.
17. The method according to claim 15 or 16, wherein the step of calculating microfracture growth is performed by repeating the following steps for all existing and new microfractures: a. checking if the location of the microfracture lies within a stress
shadow of a macrofracture of the same fracture set; b. If the location lies within such a stress shadow:
i. deactivating the microfracture; c. checking if the microfracture is active; d. If the microfracture is active:
i. increment the radius of the microfracture, based on dynamic fracture data (such as mean fracture driving stress and growth factor) for the relevant fracture set and timestep; ii. checking if the microfracture has grown into a macrofracture; iii. if the microfracture has grown into a macrofracture:
1. generating a new macrofracture of zero length at the location of the microfracture;
2. deactivating the microfracture.
18. The method according to any of claims 15-17, wherein the step of adding new macrofractures directly is performed by: a. calculating how frequently macrofractures will nucleate from growing microfractures, based on data from the implicit fracture model; b. setting current time to start of the timestep; c. repeating the following steps until the end of the timestep has been reached: i. picking a random location for a new macrofracture;
ii. checking if the location for the new macrofracture lies within a macrofracture exclusion zone;
iii. if the location does not lie within such an exclusion zone:
1. adding a new macrofracture of zero length to the
explicit discrete fracture network at the picked location, setting the status of the macrofracture tips to active; d. incrementing the current time by an amount of time corresponding to the calculated frequency of macrofractures nucleating from growing microfractures.
19. The method according to any of claims 15-18, wherein the step of calculating macrofracture growth is performed by repeating the following steps for both macrofracture tips of all macrofractures in the fracture set: a. checking if the macrofracture tip is active; b. if the macrofracture tip is active: i. calculating the maximum propagation distance of the
macrofracture tip from the dynamic data calculated for the implicit fracture model; ii. checking if the fracture tip will interact with a stress shadow of another macrofracture of the same fracture set if growing by the maximum propagation distance; iii. if such interaction will occur:
1. reducing the maximum propagation distance to the distance, at which stress shadow interaction occurs;
2. deactivating both macrofracture tips;
3. if required, adding a relay segment to link the two macrofractures. iv. checking if the fracture tip will intersect a macrofracture of another fracture set if growing by the maximum propagation distance; v. if such intersection will occur:
1. reducing the maximum propagation distance to the distance, at which intersection occurs;
2. deactivating the propagating macrofracture tip; vi. checking if the fracture tip will cross a gridblock boundary to an adjacent gridblock if growing by the maximum propagation distance; vii. if such crossing will occur:
1. calculating the real time, at which the macrofracture tip will cross the gridblock boundary;
2. reducing the maximum propagation distance to the distance, at which the gridblock boundary is reached;
3. extending the macrofracture tip by the maximum
propagation distance;
4. deactivating the propagating macrofracture tip;
5. adding a new macrofracture of zero length at the
entering position within the adjacent gridblock, connecting one of the macrofracture tips to the deactivated macrofracture tip at the gridblock boundary and setting the status of the other macrofracture tip to active;
6. checking if fracture propagation in the adjacent
gridblock has already been calculated up to this time;
7. if this is the case:
a. extending the active macrofracture tip of the new macrofracture into the adjacent gridblock in the appropriate direction of propagation until the position corresponding to the latest calculation time of the adjacent gridblock, checking for stress shadow interaction and intersection as described above in steps ii-v; viii. if no such crossing will occur:
1. extending the macrofracture tip by the maximum
propagation distance.
20. The method according to any of claims 15-19, wherein the step of
constructing the explicit discrete fracture network comprises the steps of: a. combining all macrofractures, which are connected across gridblock boundaries to create continuous macrofractures that span across multiple gridblocks; b. calculating positions of corner points of all macrofractures to ensure accurate bevelling across gridblock boundaries and macrofracture intersections.
21. The use of the method according to any of the preceding claims to generate an accurate representation of the fracture network in the subsurface, including key geometric aspects (such as fracture densities, cumulative fracture density distribution functions, fracture porosity and fracture connectivity) for an improved simulation of fluid flow through fractured rock, that can be used to manage oil and gas reservoirs, predict groundwater flow and contaminant dispersion for pollution monitoring, nuclear waste storage, predict flow and leakage in reservoirs used for carbon capture and storage (CCS) and to predict water/fluid flow rates in geothermal systems.
22. The use of the method according to any of claims 1 -20 to facilitate the
generation of accurate, risk-weighted subsurface flow models that can replicate existing oil and gas production data (history matched models), by identifying the key mechanical properties and geological parameters controlling fracture development, and generating multiple geologically plausible fracture models.
23. The use of the method according to any of claims 1 -20 to predict the impact of natural fractures on the bulk rock mechanical properties for use in assessing mine, tunnel or overland construction stability, for hydrocarbon exploration and production or for monitoring production-related subsidence or risk of seismic activity resulting from oil and gas production, CCS or CO2 injection.
24. A computer system comprising a. a processor; b. memory operative coupled to the processor, one or more modules of which memory comprises stored processor-executable instructions arranged to instruct the computer system to perform the method according to any of claims 1-20.
25. One or more non-transitory computer-readable storage media comprising computer-executable instructions arranged to instruct a computer system to perform the method according to any of claims 1-20.
PCT/EP2019/064318 2018-06-13 2019-06-03 A method and a system for modelling and simulating a fractured geological structure WO2019238451A1 (en)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
EP18177508.1 2018-06-13
EP18177508 2018-06-13

Publications (1)

Publication Number Publication Date
WO2019238451A1 true WO2019238451A1 (en) 2019-12-19

Family

ID=62630997

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/EP2019/064318 WO2019238451A1 (en) 2018-06-13 2019-06-03 A method and a system for modelling and simulating a fractured geological structure

Country Status (1)

Country Link
WO (1) WO2019238451A1 (en)

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111913216A (en) * 2020-08-03 2020-11-10 煤炭科学技术研究院有限公司 Roadway surrounding rock stability force structure cooperative monitoring method influenced by repeated mining
CN112464584A (en) * 2020-11-09 2021-03-09 长江勘测规划设计研究有限责任公司 Method for estimating water level and flow of free surface flow
CN112711884A (en) * 2020-12-29 2021-04-27 华南理工大学 Display and hidden algorithm combined brittle material quasi-static damage simulation prediction method
CN113156080A (en) * 2021-04-08 2021-07-23 青岛海洋地质研究所 Device and method for simulating influence law of diapir action on hydrate accumulation
CN113268856A (en) * 2021-04-20 2021-08-17 南方科技大学 Probability-based crack characterization and reconstruction method, storage medium and terminal device
CN113284089A (en) * 2021-04-20 2021-08-20 深圳大学 Crack generation method based on generator, storage medium and terminal equipment
CN113343510A (en) * 2021-05-06 2021-09-03 武汉大学 Two-dimensional random fracture grid generation method
CN114718446A (en) * 2022-04-18 2022-07-08 中铁二院工程集团有限责任公司 Mountain railway tunnel drilling arrangement method and deep hole drilling method
CN114970260A (en) * 2022-05-20 2022-08-30 武汉大学 Lattice phase field method for simulating composite material damage
WO2023130074A1 (en) * 2021-12-31 2023-07-06 Schlumberger Technology Corporation Geologic modeling framework
CN117150978A (en) * 2023-11-01 2023-12-01 中国地质大学（北京） Rock mass change prediction system based on rock mass fracture seepage information
CN114647004B (en) * 2022-02-25 2024-03-15 中国地质大学（北京） Method for confirming sliding directions of underground sliding fracture at different periods

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130054207A1 (en) * 2009-06-05 2013-02-28 Schlumberger Technology Corporation Fracture Network Characterization Method
US20170316128A1 (en) * 2016-04-29 2017-11-02 Hao Huang Method and system for characterizing fractures in a subsurface region

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130054207A1 (en) * 2009-06-05 2013-02-28 Schlumberger Technology Corporation Fracture Network Characterization Method
US20170316128A1 (en) * 2016-04-29 2017-11-02 Hao Huang Method and system for characterizing fractures in a subsurface region

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
BONNEAU FRANÇOIS ET AL: "A methodology for pseudo-genetic stochastic modeling of discrete fracture networks", COMPUTERS AND GEOSCIENCES, PERGAMON PRESS, OXFORD, GB, vol. 56, 24 February 2013 (2013-02-24), pages 12 - 22, XP028542830, ISSN: 0098-3004, DOI: 10.1016/J.CAGEO.2013.02.004 *
L. SOUCHE ET AL: "An End-to End-Approach to Naturally Fractured Reservoir Modeling: Workflow and Implementation", SEG/SAN ANTONIO 2007 ANNUAL MEETING, 3 September 2007 (2007-09-03), us, XP055617817, ISSN: 2214-4609, ISBN: 978-90-73781-53-5, DOI: 10.3997/2214-4609.20146718 *

Cited By (20)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111913216B (en) * 2020-08-03 2022-09-20 煤炭科学技术研究院有限公司 Roadway surrounding rock stability force structure cooperative monitoring method influenced by repeated mining
CN111913216A (en) * 2020-08-03 2020-11-10 煤炭科学技术研究院有限公司 Roadway surrounding rock stability force structure cooperative monitoring method influenced by repeated mining
CN112464584A (en) * 2020-11-09 2021-03-09 长江勘测规划设计研究有限责任公司 Method for estimating water level and flow of free surface flow
CN112464584B (en) * 2020-11-09 2023-03-24 长江勘测规划设计研究有限责任公司 Method for estimating water level and flow of free surface flow
CN112711884B (en) * 2020-12-29 2024-01-23 华南理工大学 Simulation prediction method for quasi-static damage of brittle material by combination of implicit algorithm
CN112711884A (en) * 2020-12-29 2021-04-27 华南理工大学 Display and hidden algorithm combined brittle material quasi-static damage simulation prediction method
CN113156080A (en) * 2021-04-08 2021-07-23 青岛海洋地质研究所 Device and method for simulating influence law of diapir action on hydrate accumulation
CN113156080B (en) * 2021-04-08 2022-03-25 青岛海洋地质研究所 Device and method for simulating influence law of diapir action on hydrate accumulation
CN113268856B (en) * 2021-04-20 2022-05-20 南方科技大学 Probability-based crack characterization and reconstruction method, storage medium and terminal device
CN113284089A (en) * 2021-04-20 2021-08-20 深圳大学 Crack generation method based on generator, storage medium and terminal equipment
CN113268856A (en) * 2021-04-20 2021-08-17 南方科技大学 Probability-based crack characterization and reconstruction method, storage medium and terminal device
CN113343510B (en) * 2021-05-06 2022-05-17 武汉大学 Two-dimensional random fracture grid generation method
CN113343510A (en) * 2021-05-06 2021-09-03 武汉大学 Two-dimensional random fracture grid generation method
WO2023130074A1 (en) * 2021-12-31 2023-07-06 Schlumberger Technology Corporation Geologic modeling framework
CN114647004B (en) * 2022-02-25 2024-03-15 中国地质大学（北京） Method for confirming sliding directions of underground sliding fracture at different periods
CN114718446A (en) * 2022-04-18 2022-07-08 中铁二院工程集团有限责任公司 Mountain railway tunnel drilling arrangement method and deep hole drilling method
CN114970260A (en) * 2022-05-20 2022-08-30 武汉大学 Lattice phase field method for simulating composite material damage
CN114970260B (en) * 2022-05-20 2024-04-02 武汉大学 Lattice phase field method for simulating composite material damage
CN117150978B (en) * 2023-11-01 2024-01-05 中国地质大学（北京） Rock mass change prediction system based on rock mass fracture seepage information
CN117150978A (en) * 2023-11-01 2023-12-01 中国地质大学（北京） Rock mass change prediction system based on rock mass fracture seepage information

Similar Documents

Publication Publication Date Title
WO2019238451A1 (en) A method and a system for modelling and simulating a fractured geological structure
US10712472B2 (en) Method and system for forming and using a subsurface model in hydrocarbon operations
US10846447B2 (en) Method and system for stacking fracture prediction
US10572611B2 (en) Method and system for characterizing fractures in a subsurface region
RU2669948C2 (en) Multistage oil field design optimisation under uncertainty
US8190414B2 (en) Modeling of hydrocarbon reservoirs containing subsurface features
US9068448B2 (en) System and method for predicting fluid flow characteristics within fractured subsurface reservoirs
EP2616979B1 (en) Production estimation in subterranean formations
EP1825303B1 (en) Method system and program storage device for optimization of valve settings in instrumented wells using adjoint gradient technology and reservoir simulation
US10810331B2 (en) System for predicting induced seismicity potential resulting from injection of fluids in naturally fractured reservoirs
AU2013399052B2 (en) Reservoir simulator, method and computer program product
US11789170B2 (en) Induced seismicity
CN110632657B (en) Mudstone smearing type fault sealing analysis method and device
US11434759B2 (en) Optimization of discrete fracture network (DFN) using streamlines and machine learning
Feng et al. Development characteristics and quantitative prediction of multiperiod fractures in superdeep thrust-fold belt
Jang et al. The oil production performance analysis using discrete fracture network model with simulated annealing inverse method
Welch et al. DFN Generator v2. 0: A new tool to model the growth of large-scale natural fracture networks using fundamental geomechanics
Ogbechie Fracture modeling and flow behavior in shale gas reservoirs using discrete fracture networks
CN113031056B (en) Fault closure analysis method and device under construction constraint
Cruz et al. Hydraulic fracture propagation in a vertically and laterally heterogeneous stress media in the Permian Basin
Pankaj et al. Hydraulic Fracture Calibration for Unconventional Reservoirs: A New Methodology for Predictive Modelling
Tran* XD et al. Practical Reservoir Simulation for Small Development Teams-Customizing Stimulation Designs by Landing Zone in the Midland Basin
Agarwal et al. Reservoir Characterization of the Ekofisk Field: A Giant, Fractured Chalk Reservoir in the Norwegian North Sea-Phase 1, Reservoir Description
Liu et al. Field Scale Stochastic Modeling of Fracture Networks: Combining Pattern Statistics with geomechanical criteria for fracture growth
Wang Numerical Simulation of Complex Hydraulic Fracture Development by Coupling Geo-mechanical and Reservoir Simulator

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19729682

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 19729682

Country of ref document: EP

Kind code of ref document: A1