SYSTEM AND METHOD FOR PERFORMING WELLBORE FRACTURE OPERATIONS

Abstract

A method of performing an oilfield operation about a wellbore penetrating a subterranean formation are provided. The method involves performing a fracture operation by generating fractures about the wellbore. The fractures define a hydraulic fracture network (HFN) about the wellbore. The method also involves generating a discrete fracture network (DFN) about the wellbore by extrapolating fracture data from the HFN. The DFN includes fracture branches with intersections therebetween and matrix blocks. The method also involves generating a depth of drainage through the DFN, defining production parameter(s), and performing a production operation to produce fluids from the subterranean formation based on the depth of drainage and the production parameter(s). The production operation may involve generating a flow rate through the DFN, generating a pressure profile of the DFN for an initial time based on the flow rate, and generating a production rate based on the pressure profile.

Description

BACKGROUND

The present disclosure relates generally to methods and systems for performing wellsite operations. More particularly, this disclosure is directed to methods and systems for performing fracture operations, such as investigating subterranean formations and characterizing hydraulic fracture networks in a subterranean formation.

In order to facilitate the recovery of hydrocarbons from oil and gas wells, the subterranean formations surrounding such wells can be hydraulically fractured. Hydraulic fracturing may be used to create cracks in subsurface formations to allow oil or gas to move toward the well. A formation is fractured by introducing a specially engineered fluid (referred to as “fracturing fluid” or “fracturing slurry” herein) at high pressure and high flow rates into the formation through one or more wellbore. Hydraulic fractures may extend away from the wellbore hundreds of feet in two opposing directions according to the natural stresses within the formation. Under certain circumstances, they may form a complex fracture network.

The fracturing fluids may be loaded with proppants, which are sized particles that may be mixed with the fracturing fluid to help provide an efficient conduit for production of hydrocarbons from the formation/reservoir to the wellbore. Proppant may comprise naturally occurring sand grains or gravel, man-made or specially engineered proppants, e.g. fibers, resin-coated sand, or high-strength ceramic materials, e.g. sintered bauxite. The proppant collects heterogeneously or homogenously inside the fracture to “prop” open the new cracks or pores in the formation. The proppant creates planes of permeable conduits through which production fluids can flow to the wellbore. The fracturing fluids are preferably of high viscosity, and therefore capable of carrying effective volumes of proppant material.

The fracturing fluid may be realized by a viscous fluid, sometimes referred to as a “pad” that is injected into the treatment well at a rate and pressure sufficient to initiate and propagate a fracture in hydrocarbon formation. Injection of the “pad” is continued until a fracture of sufficient geometry is obtained to permit placement of the proppant particles. After the injection of the “pad,” the fracturing fluid may consist of a fracturing fluid and proppant material. The fracturing fluid may be a gel, oil based, water based, brine, acid, emulsion, foam or any other similar fluid. The fracturing fluid can contain several additives, viscosity builders, drag reducers, fluid-loss additives, corrosion inhibitors and the like. In order to keep the proppant suspended in the fracturing fluid until such time as all intervals of the formation have been fractured as desired, the proppant may have a density close to the density of the fracturing fluid utilized.

Proppants may be comprised of any of the various commercially available fused materials such as silica or oxides. These fused materials can comprise any of the various commercially available glasses or high-strength ceramic products. Following the placement of the proppant, the well may be shut-in for a time sufficient to permit the pressure to bleed off into the formation. This causes the fracture to close and exert a closure stress on the propping agent particles. The shut-in period may vary from a few minutes to several days.

Current hydraulic fracture monitoring methods and systems may map where the fractures occur and the extent of the fractures. Some methods and systems of microseismic monitoring may process seismic event locations by mapping seismic arrival times and polarization information into three-dimensional space through the use of modeled travel times and/or ray paths. These methods and systems can be used to infer hydraulic fracture propagation over time.

Conventional hydraulic fracture models may also assume a bi-wing type induced fracture. These bi-wing fractures may be short in representing the complex nature of induced fractures in some unconventional reservoirs with preexisting natural fractures. Published models may map the complex geometry of discrete hydraulic fractures based on monitoring microseismic event distribution.

In some cases, models may not be constrained by accounting for either the amount of pumped fluid or mechanical interactions between fractures and injected fluid and among the fractures. Some of the constrained models may provide a fundamental understanding of involved mechanisms, but may be complex in mathematical description and/or require computer processing resources and time in order to provide accurate simulations of hydraulic fracture propagation.

Unconventional formations, such as shales are being developed as sources of hydrocarbon production. Once considered only as source rocks and seals, shale formations are now considered as tight-porosity and low-permeability unconventional reservoirs. Patterns of hydraulic fractures created by the fracturing stimulation may be complex and may form a fracture network as indicated by the distribution of associated microseismic events. Complex hydraulic fracture networks have been developed to represent the created hydraulic fractures. Examples of fracture models are provided in U.S. Pat. Nos. 6,101,447, 7,363,162, 7,788,074, 20080133186, 20100138196, and 20100250215.

Hydraulic fracturing of shale formation may be used to stimulate and produce from the reservoir. Production simulation has been developed to estimate production from reservoirs. Various production simulation techniques have been used with conventional reservoirs. Examples of production simulation are provided in Warren et al., “The Behavior of Naturally Fractured Reservoirs, Soc.Pet.Eng.J., Vol. 3(3): pp. 245-255 (1963) (hereafter “Warren & Root”); Basquet et al., “Gas Flow Simulation in Discrete Fracture Network Models”. Paper SPE 79708 presented at the SPE Reservoir Simulation Symposium, Houston, Tex., 3-5 Feb. 2003 (hereafter “Basquet”); Gong et al., “Detailed Modeling of the Complex Fracture Network of Shale Gas Reservoirs”, SPE paper 142705 presented at the SPE Middle East Unconventional Gas Conference and Exhibition held in Muscat, Oman, 31 Jan. 2011 (hereafter “Gong”); Cinco-Ley et al., “Pressure Transient Analysis for Naturally Fractured Reservoirs”, SPE paper 11026 presented at the Annual Fall Technical Conference and Exhibition held in New Orleans, La., Sep. 26, 1982 (hereafter “Cinco-Ley”); Xu et al., “Quick Estimate of Initial Production from Stimulated Reservoirs with Complex Hydraulic Fracture Network”, Paper SPE 146753 presented at the SPE Annual Technical Conference and Exhibition held in Denver, Colo., USA, 30 Oct. 2011 (hereafter “Xu 2011”); and C. E. Cohen et al. “Production Forecast After Hydraulic Fracturing in Naturally Fractured Reservoir: Coupling a Complex Fracturing Simulator and a Semi-Analytical Production Model”, Paper (SPE 152541) presented at the SPE Hydraulic Fracturing Technology Conference and Exhibition held in The Woodlands, Tex., USA, 8 Feb. 2012, the entire contents of which are hereby incorporated by reference. However, the reservoirs may be unconventional and/or have natural fractures, such as those with shale.

SUMMARY

In at least one aspect the present disclosure relates to a method of performing a production operation about a wellbore penetrating a subterranean formation. The subterranean formation has a plurality of fractures thereabout. The method involves generating a flow rate through a discrete fracture network defined by the plurality of fractures in the subterranean formation. The discrete fracture network includes a plurality of fracture branches with intersections therebetween and a plurality of matrix blocks. The method further involves generating a pressure profile of the discrete fracture network for an initial time based on the flow rate, and generating a production rate based on the pressure profile.

In another aspect, the disclosure relates to a method of performing an oilfield operation about a wellbore penetrating a subterranean formation. The method involves performing a fracture operation comprising generating fractures about the wellbore. The fractures define a hydraulic fracture network about the wellbore. The method also involves generating a discrete fracture network about the wellbore by extrapolating fracture data from the hydraulic fracture network. The discrete fracture network includes a plurality of fracture branches with intersections therebetween and a plurality of matrix blocks. The method further involves generating a depth of drainage through the discrete fracture network, defining at least one production parameter, and performing a production operation to produce fluids from the subterranean formation based on the depth of drainage and the at least one production parameter.

Finally, in another aspect, the disclosure relates to a method of performing an oilfield operation about a wellbore penetrating a subterranean formation. The method involves stimulating the wellbore by injecting fluid into the subterranean formation such that fractures are generated about the wellbore, measuring the fractures and defining a hydraulic fracture network based on the measured fractures.

The method also involves generating a discrete fracture network about the wellbore by extrapolating fracture data from the hydraulic fracture network. The discrete fracture network includes a plurality of fracture branches with intersections therebetween and a plurality of matrix blocks. The method also involves generating a depth of drainage through the discrete fracture network, defining at least one production parameter, estimating a production rate over time based on the depth of drainage and the production parameter(s), and producing fluids from the subterranean formation based on the estimated production rate.

This summary is provided to introduce a selection of concepts that are further described below in the detailed description. This summary is not intended to identify key or essential features of the claimed subject matter, nor is it intended to be used as an aid in limiting the scope of the claimed subject matter.

BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments of the system and method for characterizing wellbore stresses are described with reference to the following figures. The same numbers are used throughout the figures to reference like features and components.

FIGS. 1.1-1.4 are schematic views illustrating various oilfield operations at a wellsite;

FIGS. 2.1-2.4 are schematic views of data collected by the operations of FIGS. 1.1-1.4;

FIG. 3 is a schematic illustration of a hydraulic fracturing site depicting a fracture operation;

FIGS. 4.1 and 4.2 are flow charts depicting methods of performing an oilfield operation and a production operation, respectively;

FIG. 5 is a schematic illustration of a production simulation of a discrete fracture network (DFN) extracted from a hydraulic fracturing simulation;

FIG. 6 is a schematic illustration of the DFN of FIG. 5 having a plurality of matrix blocks;

FIG. 7. is a schematic illustration of a an approximation of flow through a matrix block;

FIGS. 8.1-8.3 are graphs illustrating production, cumulated production and pressure, respectively, of a well;

FIG. 9 is a schematic diagram depicting coordinates of fractures of a matrix block;

FIG. 10 is a schematic diagram depicting flow rate from a matrix block to a branch of a DFN;

FIGS. 11.1 and 11.2 are graphs depicting pressure versus time over time for a highly conductive DFN;

FIG. 12 is a graph of normalized pressure and time delay over time for a highly conductive DFN;

FIG. 13 is a graph of cumulated production over time for a highly conductive DFN;

FIGS. 14.1 and 14.2 are graphs depicting pressure versus time over time for a low conductivity DFN;

FIG. 15 a graph of normalized pressure and time delay over time for a low conductivity DFN;

FIG. 16 is a graph of cumulated production over time for a low conductivity DFN;

FIG. 17 is a graph of normalized pressure and time delay over time for a low conductivity DFN using an Unconventional Production Model (UPM);

FIG. 18 is a graph of cumulated production over time for a low conductivity DFN using a UPM;

FIG. 19 is a table of graphs of pressure and time delay over time;

FIG. 20 is a graph comparing simulated production over time using a reservoir simulator and the UPM;

FIGS. 21.1 and 21.2 are schematic diagrams depicting of a DFN as depicted by a reservoir simulator and the UPM, respectively;

FIG. 22 is a graph comparing simulated production over time for different fracture conductivities using a reservoir simulator and the UPM; and

FIGS. 23.1 and 23.2 are graphs of flow rate and cumulated production, respectively, over time by a reservoir simulator, the UPM and the UPM without delay.

DETAILED DESCRIPTION

The description that follows includes exemplary systems, apparatuses, methods, and instruction sequences that embody techniques of the subject matter herein. However, it is understood that the described embodiments may be practiced without these specific details.

The present disclosure relates to techniques for performing fracture operations to estimate and/or predict production. The fracture operations involve fracture modeling that utilize elliptical and wire mesh modeling to estimate production.

FIGS. 1.1-1.4 depict various oilfield operations that may be performed at a wellsite, and FIGS. 2.1-2.4 depict various information that may be collected at the wellsite. FIGS. 1.1-1.4 depict simplified, schematic views of a representative oilfield or wellsite 100 having subsurface formation 102 containing, for example, reservoir 104 therein and depicting various oilfield operations being performed on the wellsite 100. FIG. 1.1 depicts a survey operation being performed by a survey tool, such as seismic truck 106.1, to measure properties of the subsurface formation. The survey operation may be a seismic survey operation for producing sound vibrations. In FIG. 1.1, one such sound vibration 112 generated by a source 110 reflects off a plurality of horizons 114 in an earth formation 116. The sound vibration(s) 112 may be received in by sensors, such as geophone-receivers 118, situated on the earth's surface, and the geophones 118 produce electrical output signals, referred to as data received 120 in FIG. 1.1.

In response to the received sound vibration(s) 112 representative of different parameters (such as amplitude and/or frequency) of the sound vibration(s) 112, the geophones 118 may produce electrical output signals containing data concerning the subsurface formation. The data received 120 may be provided as input data to a computer 122.1 of the seismic truck 106.1, and responsive to the input data, the computer 122.1 may generate a seismic and microseismic data output 124. The seismic data output may be stored, transmitted or further processed as desired, for example by data reduction.

FIG. 1.2 depicts a drilling operation being performed by a drilling tool 106.2 suspended by a rig 128 and advanced into the subsurface formations 102 to form a wellbore 136 or other channel. A mud pit 130 may be used to draw drilling mud into the drilling tools via flow line 132 for circulating drilling mud through the drilling tools, up the wellbore 136 and back to the surface. The drilling mud may be filtered and returned to the mud pit. A circulating system may be used for storing, controlling or filtering the flowing drilling muds. In this illustration, the drilling tools are advanced into the subsurface formations to reach reservoir 104. Each well may target one or more reservoirs. The drilling tools may be adapted for measuring downhole properties using logging while drilling tools. The logging while drilling tool may also be adapted for taking a core sample 133 as shown, or removed so that a core sample may be taken using another tool.

A surface unit 134 may be used to communicate with the drilling tools and/or offsite operations. The surface unit may communicate with the drilling tools to send commands to the drilling tools, and to receive data therefrom. The surface unit may be provided with computer facilities for receiving, storing, processing, and/or analyzing data from the operation. The surface unit may collect data generated during the drilling operation and produce data output 135 which may be stored or transmitted. Computer facilities, such as those of the surface unit, may be positioned at various locations about the wellsite and/or at remote locations.

Sensors (S), such as gauges, may be positioned about the oilfield to collect data relating to various operations as described previously. As shown, the sensor (S) may be positioned in one or more locations in the drilling tools and/or at the rig to measure drilling parameters, such as weight on bit, torque on bit, pressures, temperatures, flow rates, compositions, rotary speed and/or other parameters of the operation. Sensors (S) may also be positioned in one or more locations in the circulating system.

The data gathered by the sensors may be collected by the surface unit and/or other data collection sources for analysis or other processing. The data collected by the sensors may be used alone or in combination with other data. The data may be collected in one or more databases and/or transmitted on or offsite. All or select portions of the data may be selectively used for analyzing and/or predicting operations of the current and/or other wellbores. The data may be historical data, real time data or combinations thereof. The real time data may be used in real time, or stored for later use. The data may also be combined with historical data or other inputs for further analysis. The data may be stored in separate databases, or combined into a single database.

The collected data may be used to perform analysis, such as modeling operations. For example, the seismic data output may be used to perform geological, geophysical, and/or reservoir engineering analysis. The reservoir, wellbore, surface and/or processed data may be used to perform reservoir, wellbore, geological, and geophysical or other simulations. The data outputs from the operation may be generated directly from the sensors, or after some preprocessing or modeling. These data outputs may act as inputs for further analysis.

The data may be collected and stored at the surface unit 134. One or more surface units may be located at the wellsite, or connected remotely thereto. The surface unit may be a single unit, or a complex network of units used to perform the necessary data management functions throughout the oilfield. The surface unit may be a manual or automatic system. The surface unit 134 may be operated and/or adjusted by a user.

The surface unit may be provided with a transceiver 137 to allow communications between the surface unit and various portions of the current oilfield or other locations. The surface unit 134 may also be provided with or functionally connected to one or more controllers for actuating mechanisms at the wellsite 100. The surface unit 134 may then send command signals to the oilfield in response to data received. The surface unit 134 may receive commands via the transceiver or may itself execute commands to the controller. A processor may be provided to analyze the data (locally or remotely), make the decisions and/or actuate the controller. In this manner, operations may be selectively adjusted based on the data collected. Portions of the operation, such as controlling drilling, weight on bit, pump rates or other parameters, may be optimized based on the information. These adjustments may be made automatically based on computer protocol, and/or manually by an operator. In some cases, well plans may be adjusted to select optimum operating conditions, or to avoid problems.

FIG. 1.3 depicts a wireline operation being performed by a wireline tool 106.3 suspended by the rig 128 and into the wellbore 136 of FIG. 1.2. The wireline tool 106.3 may be adapted for deployment into a wellbore 136 for generating well logs, performing downhole tests and/or collecting samples. The wireline tool 106.3 may be used to provide another method and apparatus for performing a seismic survey operation. The wireline tool 106.3 of FIG. 1.3 may, for example, have an explosive, radioactive, electrical, or acoustic energy source 144 that sends and/or receives electrical signals to the surrounding subsurface formations 102 and fluids therein.

The wireline tool 106.3 may be operatively connected to, for example, the geophones 118 and the computer 122.1 of the seismic truck 106.1 of FIG. 1.1. The wireline tool 106.3 may also provide data to the surface unit 134. The surface unit 134 may collect data generated during the wireline operation and produce data output 124 which may be stored or transmitted. The wireline tool 106.3 may be positioned at various depths in the wellbore to provide a survey or other information relating to the subsurface formation.

Sensors (S), such as gauges, may be positioned about the wellsite 100 to collect data relating to various operations as described previously. As shown, the sensor (S) is positioned in the wireline tool 106.3 to measure downhole parameters which relate to, for example porosity, permeability, fluid composition and/or other parameters of the operation.

FIG. 1.4 depicts a production operation being performed by a production tool 106.4 deployed from a production unit or Christmas tree 129 and into the completed wellbore 136 of FIG. 1.3 for drawing fluid from the downhole reservoirs into surface facilities 142. Fluid flows from reservoir 104 through perforations in the casing (not shown) and into the production tool 106.4 in the wellbore 136 and to the surface facilities 142 via a gathering network 146.

Sensors (S), such as gauges, may be positioned about the oilfield to collect data relating to various operations as described previously. As shown, the sensor (S) may be positioned in the production tool 106.4 or associated equipment, such as the Christmas tree 129, gathering network, surface facilities and/or the production facility, to measure fluid parameters, such as fluid composition, flow rates, pressures, temperatures, and/or other parameters of the production operation.

While only simplified wellsite configurations are shown, it will be appreciated that the oilfield or wellsite 100 may cover a portion of land, sea and/or water locations that hosts one or more wellsites. Production may also include injection wells (not shown) for added recovery or for storage of hydrocarbons, carbon dioxide, or water, for example. One or more gathering facilities may be operatively connected to one or more of the wellsites for selectively collecting downhole fluids from the wellsite(s).

It should be appreciated that FIGS. 1.2-1.4 depict tools that can be used to measure not only properties of an oilfield, but also properties of non-oilfield operations, such as mines, aquifers, storage, and other subsurface facilities. Also, while certain data acquisition tools are depicted, it will be appreciated that various measurement tools (e.g., wireline, measurement while drilling (MWD), logging while drilling (LWD), core sample, etc.) capable of sensing parameters, such as seismic two-way travel time, density, resistivity, production rate, etc., of the subsurface formation and/or its geological formations may be used. Various sensors (S) may be located at various positions along the wellbore and/or the monitoring tools to collect and/or monitor the desired data. Other sources of data may also be provided from offsite locations.

The oilfield configuration of FIGS. 1.1-1.4 depict examples of a wellsite 100 and various operations usable with the techniques provided herein. Part, or all, of the oilfield may be on land, water and/or sea. Also, while a single oilfield measured at a single location is depicted, reservoir engineering may be utilized with any combination of one or more oilfields, one or more processing facilities, and one or more wellsites.

FIGS. 2.1-2.4 are graphical depictions of examples of data collected by the tools of FIGS. 1.1-1.4, respectively. FIG. 2.1 depicts a seismic trace 202 of the subsurface formation of FIG. 1.1 taken by seismic truck 106.1. The seismic trace may be used to provide data, such as a two-way response over a period of time. FIG. 2.2 depicts a core sample 133 taken by the drilling tools 106.2. The core sample may be used to provide data, such as a graph of the density, porosity, permeability or other physical property of the core sample over the length of the core. Tests for density and viscosity may be performed on the fluids in the core at varying pressures and temperatures. FIG. 2.3 depicts a well log 204 of the subsurface formation of FIG. 1.3 taken by the wireline tool 106.3. The wireline log may provide a resistivity or other measurement of the formation at various depts. FIG. 2.4 depicts a production decline curve or graph 206 of fluid flowing through the subsurface formation of FIG. 1.4 measured at the surface facilities 142. The production decline curve may provide the production rate Q as a function of time t.

The respective graphs of FIGS. 2.1, 2.3, and 2.4 depict examples of static measurements that may describe or provide information about the physical characteristics of the formation and reservoirs contained therein. These measurements may be analyzed to define properties of the formation(s), to determine the accuracy of the measurements and/or to check for errors. The plots of each of the respective measurements may be aligned and scaled for comparison and verification of the properties.

FIG. 2.4 depicts an example of a dynamic measurement of the fluid properties through the wellbore. As the fluid flows through the wellbore, measurements are taken of fluid properties, such as flow rates, pressures, composition, etc. As described below, the static and dynamic measurements may be analyzed and used to generate models of the subsurface formation to determine characteristics thereof. Similar measurements may also be used to measure changes in formation aspects over time.

Oilfield Operations

The production operations may be simulated before, during or after production is generated from a wellbore. Simulating the production from complex fractured reservoir may be performed using various techniques. Dual porosity models may be used to address differences of properties between the fracture and the rest of the reservoir (matrix). Dual porosity may consider two coarse grids connected to each other, one for the fracture network and another one for the matrix. This method may also involve averaging of properties (e.g., for the fracture network) and simplifications to model the exchange term between the two medium. This method may be used, for example, for naturally fractured reservoirs. Additional analysis may be provided for near wellbore effects of the fracture network, such as in cases with networks created by hydraulic fracturing. Dual porosity techniques are described in Warren & Root, previously incorporated by reference herein.

Another approach involves using one medium that contains both the fracture and the reservoir, and a refined numerical grid. Additional computational time may be needed for processing. Flexibility on the gridding (e.g., unstructured mesh generation) may be provided using, for example, a specialized reservoir simulator.

Yet another approach involves the use of dual porosity equations on a discrete-fracture-network (DFN). An example of DFN is provided in Basquet, previously incorporated by reference herein. Additional methods may be used to simulate the flow from the matrix to the fractures. In some cases, such as with compressive reservoir fluids (e.g., gas), the production history from each matrix block into the DFN may be considered. The matrix block may be gridded using additional unknowns in the system of equations. Examples of gridding are provided by Gong, previously incorporated by reference herein. Analytical solutions may also be provided to simulate the flow. Solutions may be derived from a Laplace transform of the continuity equation. Examples of analytical solutions are provided by Cinco-Ley and Xu 2011, previously incorporated by reference herein.

Transient fracture pressure may be considered to obtain a complex expression that may use numerical integration in time. Constant fracture pressure may also be considered, and an expression of the flow rate between the matrix and the fracture that is linear in pressure may also be obtained. This solution may be used, for example, in conductive fractures where the variations of pressure inside the DFN are negligible (e.g., at constant wellbore pressure). The present disclosure may employ one or more of the approaches to generate an analytical solution. This solution may extend over a range of fracture conductivity, in cases of hydraulic fracturing, and/or in naturally fractured reservoirs.

The present disclosure provides an analytical solution over a range of fracture conductivity in cases of hydraulic fracturing in naturally fractured reservoirs. Such simulation may apply to unconventional reservoirs, such as shale gas, although it can be applicable to other subterranean formations as well. These unconventional reservoirs have two main features: low rock permeability and a dense network of natural fractures. A stimulation approach may be provided to address potential differences in the production mode of unconventional or other reservoirs which may involve horizontal wells and large hydraulic fracturing treatment to produce. In some cases, these treatments initiate hydraulic fractures that interact with natural fractures, and may result in complex fracturing network that connect the well to the reservoir.

This disclosure discloses a methodology to simulate the production from reservoirs, such as unconventional (naturally fractured) reservoir after a complex network of hydraulic fractures has been created. The disclosed method first extrapolates the results from an unconventional fracture model (UFM) simulation and then process it with a methodology that would give to the user a forecast of the production of the well for several years, within a time limit and accuracy range. The method of the current application extends the validity of the semi-analytical model for a full range of fracture conductivities to consider in real cases. The simulator may be validated against simulations by reservoir simulations, such as ECLIPSE™ commercially available from Schlumberger Technology Corporation (see: www.slb.com), to illustrate the capabilities of the algorithm to provide accurate results for a given range of fracture conductivity.

The current disclosure also discloses a method to simulate the production from a naturally fractured reservoir that has been stimulated by hydraulic fracturing. Portions of the method may be implemented into a software program that simulates hydraulic fracture treatments. The method may first extrapolate the results from the simulation to re-create an adapted hydraulic fracture network with averaged properties between intersections of the network, and then estimate equivalent block depth in front of each fracture face. Finally, parameters may be input for the production condition and the production simulator run. The production from each matrix block in contact with a fracture uses an analytic expression that can be extended to a full range of realistic values for the problem parameters (conductivity, permeability, etc.) This is accomplished by updating at each time step and for each fracture face, the initial production time in order to accounts for delays in the production of each matrix block, and in a way to preserve mass balance of reservoir fluid in place/produced. The update is carried by a search algorithm that calculates this initial time so the mass really produced from each side of a matrix block is equal to the same mass if the current pressure conditions in the adjacent fracture would have been constant in time and had started at the updated initial time. The method may be compared with simulations by existing reservoir simulators, such as ECLIPSE™. Results over a range of fracture conductivity may be performed and cross-checked with a reservoir simulator.

FIG. 3 illustrates an exemplary operational setting for hydraulic fracturing of a subterranean formation (referred to herein as a “fracture site”) in accordance with the present disclosure. The fracture site 300 can be located on land or in a water environment and includes a treatment well 301 extending into a subterranean formation as well as a monitoring well 303 extending into the subterranean formation and offset from the treatment well 301. The monitoring well 303 includes an array of geophone receivers 305 (e.g., three-component geophones) spaced therein as shown.

During the fracturing operation, fracturing fluid is pumped from the surface 311 into the treatment 301 causing the surrounding formation in a hydrocarbon reservoir 307 to fracture and form a hydraulic fracture network 308. Such fracturing produces microseismic events 310, which emit both compressional waves (also referred to as primary waves or P-waves) and shear waves (also referred to as secondary waves or S-waves) that propagate through the earth and are recorded by the geophone receiver array 305 of the monitoring well 303.

The distance to the microseismic events 310 can be calculated by measuring the difference in arrival times between the P-waves and the S-waves. Also, hodogram analysis, which examines the particle motion of the P-waves, can be used to determine azimuth angle to the event. The depth of the event 310 is constrained by using the P- and S-wave arrival delays between receivers of the array 305. The distance, azimuth angle and depth values of such microseismic events 310 can be used to derive a geometric boundary or profile of the fracturing caused by the fracturing fluid over time, such as an elliptical boundary defined by a height h, elliptic aspect ratio e and major axis a as illustrated in FIG. 3.

The site 301 also includes a supply of fracturing fluid and pumping means (not shown) for supplying fracturing fluid under high pressure to the treatment well 301. The fracturing fluid can be stored with proppant (and possibly other special ingredients) pre-mixed therein. Alternatively, the fracturing fluid can be stored without pre-mixed proppant or other special ingredients, and the proppant (and/or other special ingredients) mixed into the fracturing fluid in a controlled manner by a process control system as described in U.S. Pat. No. 7,516,793, herein incorporated by reference in its entirety. The treatment well 301 also includes a flow sensor S as schematically depicted for measuring the pumping rate of the fracturing fluid supplied to the treatment well and a downhole pressure sensor for measuring the downhole pressure of the fracturing fluid in the treatment well 301.

A data processing system 309 is linked to the receivers of the array 305 of the monitoring well 303 and to the sensor S (e.g., flow sensor and downhole pressure sensor) of the treatment well 301. The data processing system 309 may be incorporated into and/or work with the surface unit 134. The data processing system 309 carries out the processing set forth in FIG. 4 and described herein. As will be appreciated by those skilled in the art, the data processing system 309 includes data processing functionality (e.g., one or more microprocessors, associated memory, and other hardware and/or software) to implement the disclosure as described herein.

The data processing system 309 can be realized by a workstation or other suitable data processing system located at the site 301. Alternatively, the data processing system 309 can be realized by a distributed data processing system wherein data is communicated (preferably in real time) over a communication link (typically a satellite link) to a remote location for data analysis as described herein. The data analysis can be carried out on a workstation or other suitable data processing system (such as a computer cluster or computing grid). Moreover, the data processing functionality of the present disclosure can be stored on a program storage device (e.g., one or more optical disks or other hand-holdable non-volatile storage apparatus, or a server accessible over a network) and loaded onto a suitable data processing system as needed for execution thereon as described herein.

FIG. 4.1 is a flow chart depicting a method 400.1 of performing an oilfield operation. The method involves 420 performing a fracture operation (actual or simulated), 422 generating a DFN about the wellbore, 424 generating a depth of drainage through the DFN, 426 defining at least one production parameter, and 428 performing a production operation.

FIG. 4.2 depicts a method 400.2 of performing a production operation. This production operation may be the same as the production operation 428 of FIG. 4.1. In the version of method 400.2, the production operation is simulated. As indicated in FIG. 4.2, the method 400.2 involves 421 generating a flow rate through a discrete fracture network, 423 generating a pressure profile of the discrete fracture network based on the flow rate, and 425 generating a production rate based on the pressure profile. The method may also involve 427 validating the production rates. The method may be provided with other features, and performed in any order.

The performing a fracture operation 420 involves generating fractures about the wellbore and defining a hydraulic fracture network about the wellbore. This fracture operation may be performed by actual injection of fluid as shown, for example, in FIG. 3. Hydraulic fracturing of the well may also be simulated using hydraulic fracture simulators. Simulations may involve generating a fracture network about the wellbore. Discrete fracture network techniques are provided in US Patent Application No. 20100307755. Data from the actual or simulated hydraulic fracturing may be used to generate data describing the resulting DFN.

A hydraulic fracture simulation 530 may be visually depicted by computer generated images as shown in FIG. 5. The hydraulic fracture simulation 530 includes a plurality of fractures 534 that form a hydraulic fracture network 536. Features of the fracture network 536, such as slurry 538, fluid 540 and bank 542, are depicted in the fracture network 536.

Generating a DFN 422 involves extrapolating fracture data from the hydraulic fracture network. The DFN may be generated by extrapolation of fracture data. The fracture data may be extrapolated from the hydraulic fracture simulation 530. This data may be exported automatically to form a production network visualization 532 as schematically depicted by arrow 533. FIG. 5 shows an example of a data export from the hydraulic fracture simulation 530 to the production network visualization 532. The production network visualization 532 provides an example of the creation of a simulated hydraulic fracture network from measured fracture data into an equivalent DFN network. The export may be performed to create a DFN 535 in a format that can be used by a production model.

In the example shown in FIG. 5, the DFN 535 includes branches 544 and intersections (or fracture tips) 546. These fracture branches 544 and intersections 546 extract portions of the hydraulic fracture simulation 530 that depict fluid flow through the fracture network 536. The remainder of the fractures 534 has been eliminated.

The format of the DFN 535 considers a unique averaged value for each property at each fracture branch 544. The fracture branches 544 are defined as the plane that connects two intersections 536. These intersections 536 may be a fracture intersection, or a fracture intersection and a fracture tip. The properties at each fracture branch 544 may be, for example, spatial coordinates at an extremity of the branch, the averaged conductivity, the averaged height, the averaged reservoir pressure at the branch location, and/or the averaged reservoir permeability at the branch location.

The description of the DFN 535 by the intersections 546 and branches 544 may be used by the present model to calculate the pressure at the intersections 546. This description may also use the branches 544 to both connect the intersections 546 and to calculate production from adjacent matrix blocks.

Referring back to FIG. 4, the generating 424 a depth of drainage through the DFN 535 may be carried out using matrix blocks. As shown in FIG. 6, the production network visualization 532 of FIG. 5 has been modified to a production network visualization 532 depicting a modified DFN 535′ with a matrix block 648 in front of each fracture branch 544 as shown in FIG. 6. Each of the matrix blocks has a depth 650.

The production network visualization 532′ provides an example of the generation of matrix depth 650 to be depleted on each side of all branches 544. The modified DFN 535′ may be used to automatically or manually generate a depth of drainage 650 for each matrix block 648. This may be done in such way that the total and actual volume of a given matrix block (not in contact with any reservoir boundaries) may be drained.

FIG. 7 schematically depicts flow of fluid through a matrix block. This figure illustrates the definition of equivalent block length, and the calculation of equivalent block depth. In the example shown, for a square matrix block 648.1 surrounded by four fracture branches 544 of equal length, an assumption may be made that each quarter 752 of the matrix block 648 can be depleted by the fracture branch 544 it contacts. A volume 754.1 of the matrix block 648 that is depleted and an equivalent block depth 755 to be depleted depicted.

Assuming a linear flow from the matrix block 648.1 to the fracture branch 544 (as will be described in further detail herein), it may also be assumed that this quarter 752 of the matrix block 648 has the length L of the fracture branch 544. Therefore, the depth of this “quarter” 752 of the matrix block 648.2 has to be equal to one-fourth of the block length L (or L/4) for the total volume to be depleted to be the same. As indicated by the arrow 733, using linear flow approximation, an equivalent block depth L/4 for a volume 754.2 of the matrix block 648.2 to be depleted may be determined. More complicated block shapes may be used, but may involve techniques that are more complex.

Referring again to FIG. 4, the defining 426 one or more production parameters may be performed by obtaining user input. The user may define one or more production parameters for consideration in the simulation. The user may select such production parameters based on some criteria or as desired. Examples of production parameters that may be selected include bottom hole pressure (BHP), reservoir fluid viscosity at reservoir conditions, reservoir fluid compressibility at reservoir conditions, and duration over which production is to be simulated, among others.

The performing a production operation 428 involves producing fluids from the subterranean formation based on the depth of drainage and the at least one production parameter. The production operation may be actual or simulated. Actual production operations involve producing fluids to the surface as shown in FIG. 1.4. Simulated production may be performed using production simulations. Visualization of the production results may also be provided. Such visualization may allow a user to visualize production decline and cumulated production, but also dynamics of pressure fields in the fracture network and the matrix blocks. FIGS. 8.1-8.3 provide examples of visualization of the production data in time (e.g., 140 days).

FIG. 8.1 is a graph 800.1 depicting production rate 856.1. The graph 800.1 plots production per day (Mscf/d) (y-axis) versus time t in days (x-axis). FIG. 8.2 is a graph 800.2 depicting cumulated production 856.2. The graph 800.2 depicts cumulated production P (MMscf) (y-axis) versus time t (x-axis). FIG. 8.3 is a three dimensional graph 800.3 depicting reservoir pressure (z-axis) versus distance x (m) (x-axis) and distance y (m) (y-axis), and pressure in the fracture network 858 and in the matrix blocks 848. These and other depictions may be provided. The production operation may be adjusted based on the production estimates.

Production Operations

The production operation (428 and/or 400.2) will be described in three parts. First, equations used in the analysis and their analytical solutions are presented. Second, the effect of conductivity on the model is provided, together with an example involving a single fracture branch for both high and low conductivities. Third, validity of the model and ways to address issues, such as conductivity, are provided.

1. Analytical Solution

Production rate may be determined using governing equations and analytical solutions. The continuity equation for compressible fluid in porous media is applied to both the matrix and the fractures. Inside the fracture network, the continuity equation can be re-written as follows:

$\begin{array}{cc}\frac{\partial }{\partial x_f}\left(\rho Q_f\left(x_f,t\right)\right)=-\rho Q_{\mathrm{mf}}\left(x_f,t\right)& \left(1\right)\end{array}$

Q_mf is the flow rate from the matrix to the reservoir, Q_f is the flow rate inside the fracture, ρ is the fluid density and x_f is the axis along the fracture. It is assumed that the fracture permeability (conductivity divided by the width) is so large that the transient term of the continuity equation may be neglected over the time scale considered for production simulation (from days to years). Darcy flow inside the fracture network may also be assumed.

$\begin{array}{cc}\rho Q_{\mathrm{mf}}\left(x_f,t\right)=-\frac{M}{2\mathrm{RT}}\frac{\partial }{\partial x_f}\left(-\mathrm{HC}\frac{\partial m\left(P_f\right)}{\partial x_f}\right)& \left(2\right)\end{array}$

P_f is the pressure inside the fracture, C the conductivity, T the temperature. The function m is the pseudo pressure. See Al-Hussainy et al., “The Flow of Real Gases Through Porous Media”, Journal of Petroleum Technology, 1966, pp. 624-36.

$\begin{array}{cc}m\left(P_f\right)=2{\int }_{P_{\mathrm{LB}}}^{P_f}\frac{P}{\mu \left(P\right)Z\left(P\right)}P& \left(3\right)\end{array}$

Inside the matrix, the continuity equation for a compressible fluid takes the following form.

$\begin{array}{cc}\frac{\partial }{\partial x}\left(\frac{k_mP_m}{\mu \left(P_m\right)Z\left(P_m\right)}\frac{\partial P_m}{\partial x}\right)=\frac{{\varphi }_mc_t\left(P_m\right)}{Z\left(P_m\right)}P_m\frac{\partial P_m}{\partial t}& \left(4\right)\end{array}$

P_m is the pressure inside the matrix, k_m is the matrix permeability, c_t is the fluid compressibility, μ is the viscosity, Z is the volume factor and φ_m is the porosity of the matrix. For simplification, Eq. 4 can be re-written as follows:

$\begin{array}{cc}\frac{{\partial }^2m}{\partial x_m^2}=\frac{1}{a}\frac{\partial m}{\partial t}& \left(5\right)\end{array}$

a is defined in Eq. 6.

$\begin{array}{cc}a=\frac{k_m}{{\varphi }_m\mu \left(P_m\right)c_t\left(P_m\right)}& \left(6\right)\end{array}$

To calculate Q_mf, Eq. (5) may be solved, with x_m the coordinate along an axis orthogonal to the fracture 964 and its coordinate x_f. FIG. 9 provides an illustration of the coordinates in the fracture 964 and matrix block 648.

The solution of Eq. 5 may be found by using a Laplace transform, such as explained in Jeannot, Yves. “Thansfert Thermique”, Textbook, Ecole des Mines de Nancy, 2009. http://www.thermique55.com/principal/thermique.pdf; and Bello, R. O., “Rate Transient Analysis in Shale Gas Reservoirs with Transient Linear Behavior”, PhD Thesis, 2009. A detailed list of equations, implementation, algorithm and variables for the presented method are provided herein.

The pressure profile in the matrix may be determined for constant fracture pressure. The first assumption of the model is that the gas behavior can be described by the following real gas equation:

$\begin{array}{cc}\rho =\frac{M}{\mathrm{RT}}\left(\frac{P_m}{Z\left(P_m\right)}\right)\mathrm{inside}\mathrm{the}\mathrm{matrix}& \left(7\right)\\ \rho =\frac{M}{\mathrm{RT}}\left(\frac{P_f}{Z\left(P_f\right)}\right)\mathrm{inside}\mathrm{the}\mathrm{fracture}\mathrm{network}& \left(8\right)\end{array}$

The basis equation for the linear gas flow inside the matrix block is

$\begin{array}{cc}\frac{\partial }{\partial x_m}\left(\frac{P_m}{\mu \left(P_m\right)Z\left(P_m\right)}\frac{\partial P_m}{\partial x_m}\right)=\frac{P{\varphi }_mc_t\left(P_m\right)}{Z\left(P_m\right)k_m}\frac{\partial P_m}{\partial t}\mathrm{with}& \left(9\right)\\ c_t\left(P_m\right)=\frac{1}{P_m}-\frac{1}{Z\left(P_m\right)}\frac{\partial Z\left(P_m\right)}{\partial P_m}& \left(10\right)\end{array}$

The following definition of a pseudo pressure will simplify the resolution previous equation

$\begin{array}{cc}m\left(P_m\right)=2{\int }_{P_{\mathrm{LB}}}^{P_m}\frac{P}{\mu \left(P\right)Z\left(P\right)}P& \left(11\right)\end{array}$

Equation (11) then becomes

$\begin{array}{cc}\frac{\partial m_{*}^2}{\partial x_m^2}=\frac{1}{a}\frac{\partial m_{*}}{\partial t}\mathrm{with}& \left(12\right)\\ a=\frac{k_m}{{\varphi }_m\mu \left(P_m\right)c_t\left(P_m\right)}& \left(13\right)\\ m_{*}=m\left(P_m\right)-m\left(P_{\mathrm{m\_}0}\right)=m\left(P_m\right)-m_{\mathrm{m\_}0}& \left(14\right)\end{array}$

and boundary conditions:

$\begin{array}{cc}{m}_{*}\left(x_m=L,t\right)=m_{*}\left(x_m=-L,t\right)=m_f-m_{\mathrm{m\_o}}& \left(15\right)\\ \frac{\partial m_{*}}{\partial x_m}\left(0,t\right)=0& \left(16\right)\\ m_{*}\left(x_m,t_0\right)=0& \left(17\right)\end{array}$

The Laplace transform of equation (15) gives

$\begin{array}{cc}\frac{\partial {\theta }^2}{\partial x_m^2}-q^2\theta =0\mathrm{with}& \left(18\right)\\ q^2=\frac{s}{a}& \left(19\right)\end{array}$

For the solution as the form

θ(x_m,s)=A cos h(qx_m)+B sin h(qx_m)  (20)

the Laplace transform of equation (16) is

$\begin{array}{cc}\frac{\partial \theta }{\partial x_m}\left(0,s\right)=0& \left(21\right)\end{array}$

that gives B=0, so

θ(x_m,s)=A cos h(qx_m)  (22)

For now it is assumed that the pressure in the fracture network is almost constant.

m _f(x_f,t)=m_f(x_f)

Therefore, the Laplace transform of equation (15) gives

$\begin{array}{cc}\theta \left(x_m=L,s\right)=\theta \left(x_m=-L,s\right)=\frac{m_f-m_{m_o}}{s}& \left(23\right)\end{array}$

that gives

$\begin{array}{cc}A=\frac{m_f-m_{\mathrm{m\_o}}}{s\cosh \left(\mathrm{qL}\right)}\mathrm{and}& \left(24\right)\\ \theta \left(x_m,s\right)=\frac{m_f-m_{\mathrm{m\_o}}}{s\cosh \left(\mathrm{qL}\right)}\cosh \left({\mathrm{qx}}_m\right)=\frac{m_f-m_{\mathrm{m\_o}}}{s}\frac{{}^{\mathrm{qx}}+{}^{-\mathrm{qx}}}{e^{\mathrm{qL}}\left(1+{}^{-2\mathrm{qx}}\right)}& \left(25\right)\end{array}$

By using a Taylor series development of

$\frac{1}{1+{}^{-2{\mathrm{qx}}_m}},$

the following can be provided:

$\begin{array}{cc}\theta \left(x_m,s\right)=\frac{m_f-m_{\mathrm{m\_o}}}{s}\left({}^{-q\left(L-x_m\right)}+{}^{-q\left(L+x_m\right)}\right)\sum _{n=0}^{\infty }{\left(-1\right)}^n{}^{-2\mathrm{nqL}}& \left(26\right)\end{array}$

The inverse Laplace Transform gives

$\begin{array}{cc}m\left(P_m\right)-m_{\mathrm{m\_}0}=\left(m_f-m_{\mathrm{m\_o}}\right)\sum _{n=0}^{\infty }{\left(-1\right)}^n\left(\mathrm{erfc}\left[\frac{\left(2n+1\right)L-x_m}{2\sqrt{a\left(t-t_0\right)}}\right]+\mathrm{erfc}\left[\frac{\left(2n+1\right)L+x_m}{2\sqrt{a\left(t-t_0\right)}}\right]\right)& \left(27\right)\end{array}$

Flow rate from the matrix to fracture with constant fracture pressure may then be determined. Flow rate from the matrix to the fracture is given by Darcy's law:

$\begin{array}{cc}\rho Q_{\mathrm{mf}}\left(x_f,t\right)=-\frac{{\mathrm{MHk}}_m}{2\mathrm{RT}}\sum _{k=1}^2\frac{\partial m\left(P_m\right)}{\partial x_m}{|}_{x_m=-L_k,L=L_k}& \left(28\right)\end{array}$

L_k corresponds to the maximum length of drainage on the side k of the fracture.

$\begin{array}{cc}\frac{\partial m\left(P_m\right)}{\partial x_m}=\frac{\left(m_f-m_{\mathrm{m\_o}}\right)}{\sqrt{\pi a\left(t-t_0\right)}}\sum _{n=0}^{\infty }{\left(-1\right)}^n\left({}^{-{\left(\frac{\left(2n+1\right)L-x_m}{2\sqrt{a\left(t-t_0\right)}}\right)}^2}-{}^{-{\left(\frac{\left(2n+1\right)L+x_m}{2\sqrt{a\left(t-t_0\right)}}\right)}^2}\right)& \left(29\right)\end{array}$

which gives

$\begin{array}{cc}\rho Q_{\mathrm{mf}}\left(x_f,t\right)=\left(m_f-m_{\mathrm{m\_o}}\right)\frac{{\mathrm{MHk}}_m}{2\mathrm{RT}\sqrt{\pi a}}\sum _{k=1}^2\frac{1}{\sqrt{\left(t-t_{0,k}\right)}}\sum _{n=0}^{\infty }{\left(-1\right)}^n\left({}^{-{\left(\frac{\left(n+1\right)L_k}{\sqrt{a\left(t-t_{0,k}\right)}}\right)}^2}-{}^{-{\left(\frac{nL_k}{\sqrt{a\left(t-t_{0,k}\right)}}\right)}^2}\right)& \left(30\right)\end{array}$

Flow inside the fracture branch between intersection i and j may also be determined. The flow inside the fracture network is described by the following equation:

$\begin{array}{cc}\frac{\partial }{\partial x_f}\left(\rho Q_f\left(x_f,t\right)\right)=-\rho Q_{\mathrm{mf}}\left(x_f,t\right)& \left(31\right)\end{array}$

with Q_mf the flow rate (m³/s) from the matrix to the fracture and Q_f the flux (m²/s) from the fracture. Under the assumption that the gas behavior can be described by the following real gas equation, this equation becomes:

$\begin{array}{cc}\frac{M}{2\mathrm{RT}}\frac{\partial }{\partial x_f}\left(-\mathrm{HC}\frac{\partial m\left(P_f\right)}{\partial x_f}\right)=-\rho Q_{\mathrm{mf}}\left(x_f,t\right)& \left(32\right)\end{array}$

with the following boundary conditions:

m(P_f(x_f=0))−m_m—o=m_f(x_f=0)−m_—o=m_f,i−m_—o=m_f*,i

m(P_f(x_f=L_f))−m_—o=m_f(x_f=L_f)−m_—o=m_f,j−m_—o=m_f*,j  (33)

with L_f the length of the fracture between two intersections. Using the following:

m _f* =m _f(x_f)−m_—o  (34)

and introducing equation (30) into equation (32), the following is obtained:

$\begin{array}{cc}\frac{\partial }{\partial x_f}\left(\frac{\partial m_{f*}}{\partial x_f}\right)-{\gamma }^2m_{f*}=0\mathrm{with}& \left(35\right)\\ {\gamma }^2\frac{k_m}{C\sqrt{\pi a}}\sum _{k=1}^2\frac{1}{\sqrt{\left(t-t_{0,k}\right)}}\sum _{n=0}^{\infty }{\left(-1\right)}^n\left({}^{-{\left(\frac{{\mathrm{nL}}_k}{\sqrt{a\left(t-t_{0,k}\right)}}\right)}^2}-{}^{-{\left(\frac{\left(n+1\right)L_k}{\sqrt{a\left(t-t_{0,k}\right)}}\right)}^2}\right)& \left(36\right)\end{array}$

Solution to equation (35) has the form

m _f* =Ae ^{γx f} +Be ^{−γx f}   (37)

and equations (34) gives

$\begin{array}{cc}{m}_f\left(x_f\right)=m_{\mathrm{m\_o}}+\left(\frac{m_{f*,i}{}^{-\gamma L_f}-m_{f*,j}}{{}^{-\gamma L_f}-{}^{\gamma L_f}}\right){}^{\gamma x_f}+\left(\frac{m_{f*,j}-m_{f*,i}{}^{\gamma L_f}}{{}^{-\gamma L_f}-{}^{\gamma L_f}}\right){}^{-\gamma x_f}& \left(38\right)\end{array}$

Flow rate may also be determined at intersection, for example, i from branch i,j.

$\begin{array}{cc}\rho Q_{f,i,\mathrm{ij}}=-\frac{{\mathrm{MHC}}_{i,j}}{2\mathrm{RT}}\frac{\partial m_f}{\partial x_f}{|}_{x_f=0}& \left(39\right)\end{array}$

When introducing equation (48) in equation (39)

$\begin{array}{cc}\rho Q_{f,i,\mathrm{ij}}=-\frac{\gamma {\mathrm{MHC}}_{i,j}}{2\mathrm{RT}}\left[\frac{m_{f*,i}\left({}^{-\gamma L_f}+{}^{\gamma L_f}\right)-2m_{f*,j}}{{}^{-\gamma L_f}-{}^{\gamma L_f}}\right]& \left(40\right)\end{array}$

if this is an element from the tubing, the equation becomes

$\begin{array}{cc}\rho Q_{f,i,\mathrm{ij}}=-\frac{{\mathrm{MHC}}_{i,j}}{2\mathrm{RT}}\left[\frac{m_{f*,i}-m_{f*,j}}{L_f}\right]& \left(41\right)\end{array}$

Mass balance may also be determined at an intersection between fractures

$\begin{array}{cc}\sum _{j=1}^{\mathrm{N\_ij}}\rho Q_{f,i,\mathrm{ij}}=0& \left(42\right)\end{array}$

with N_ij the number of branches reaching this intersection.

$\begin{array}{cc}\sum _{j=1}^{\mathrm{N\_ij}}\rho Q_{f,i,\mathrm{ij}}=\sum _{j=1}^{\mathrm{N\_ij}}-\frac{{\gamma }_{\mathrm{ij}}H_{\mathrm{ij}}C_{\mathrm{ij}}M}{2\mathrm{RT}}\left[\frac{\begin{array}{c}\left(m_{f,i}-m_{m_o,\mathrm{ij}}\right)\left({}^{-\gamma L_{f,\mathrm{ij}}}+{}^{\gamma L_{f,\mathrm{ij}}}\right)-\\ 2\left(m_{f,j}-m_{m_o,\mathrm{ij}}\right)\end{array}}{{}^{-\gamma L_{f,\mathrm{ij}}}-{}^{\gamma L_{f,\mathrm{ij}}}}\right]=0& \left(43\right)\end{array}$

This can be rearranged as follows

$\begin{array}{cc}{m}_{f,i}\sum _{j=1}^{\mathrm{N\_ij}}\frac{{\gamma }_{\mathrm{ij}}H_{\mathrm{ij}}C_{\mathrm{ij}}\left({}^{-\gamma L_{f,\mathrm{ij}}}+{}^{\gamma L_{f,\mathrm{ij}}}\right)}{{}^{-\gamma L_{f,\mathrm{ij}}}-{}^{\gamma L_{f,\mathrm{ij}}}}=\sum _{j=1}^{\mathrm{N\_ij}}m_{f,j}\frac{{\gamma }_{\mathrm{ij}}H_{\mathrm{ij}}C_{\mathrm{ij}}}{{}^{-\gamma L_{f,\mathrm{ij}}}-{}^{\gamma L_{f,\mathrm{ij}}}}=\sum _{j=1}^{\mathrm{N\_ij}}m_{m_o,\mathrm{ij}}\frac{{\gamma }_{\mathrm{ij}}H_{\mathrm{ij}}C_{\mathrm{ij}}\left({}^{-\gamma L_{f,\mathrm{ij}}}+{}^{\gamma L_{f,\mathrm{ij}}}-2\right)}{{}^{-\gamma L_{f,\mathrm{ij}}}-{}^{\gamma L_{f,\mathrm{ij}}}}& \left(44\right)\end{array}$

The time function t_o,k(t) of fluid flow through the matrix block may also be updated. Objective function F may be defined as the difference between the real mass produced so far from each face of each matrix block, and the mass that would have been produced if the current pressure field inside the DFN had been constant and the initial time considered in the analytical solution had been constant and equal to t_o,k(t).

$\begin{array}{cc}{F}_k\left(t,t_{0,k}\left(t\right)\right)={\int }_0^t\rho Q_{\mathrm{tot},\mathrm{ij},k}\left(\tau ,P\left(\tau \right),t_{0,k}\left(\tau \right)\right)\tau -{\int }_{t_{0,k}\left(t\right)}^t\rho Q_{\mathrm{tot},\mathrm{ij},k}\left(\tau ,P\left(t\right),t_{0,k}\left(t\right)\right)\tau \mathrm{or}F_k\left(t,t_{0,k}\left(t\right)\right)=M_{\mathrm{tot},k}\left(t\right)-{\int }_{t_{0,k}\left(t\right)}^t{\int }_0^{L_f}\rho Q_{\mathrm{mf},k}\left(\tau ,P_f\left(x_f,t\right),t_{0,k}\left(t\right)\right)\tau x_f& \left(45\right)\end{array}$

Initial time t_o,k(t) may be calculated by finding the value of t_o,k(t) that would give F_o,k(t)=0. Thus, the total mass produced from face k of the branch is equal to the mass that would have been produced so far by the same branch under certain conditions, such as if the current pressure condition in the fracture would have been the same and constant from initial production time t_o,k(t), and/or if there would have been no production from this face before the initial time t_o,k(t).

This search for the value of t_o,k(t) is done with the iterative algorithm of Newton:

$\begin{array}{cc}{t}_{0,k}^{n+1}\left(t\right)=t_{0,k}^n\left(t\right)-\frac{F_k^n\left(t,t_{0,k}^n\left(t\right)\right)}{F_k^{\prime ,n}\left(t,t_{0,k}^n\left(t\right)\right)}\mathrm{with}F_k^{\prime ,n}\left(t,t_{0,k}^n\left(t\right)\right)=\frac{\partial F_k^n\left(t,t_{0,k}^n\left(t\right)\right)}{\partial t_{0,k}^n\left(t\right)}& \left(46\right)\end{array}$

The derivative

$\frac{\partial F_k^n\left(t,t_{0,k}^n\left(t\right)\right)}{\partial t_{0,k}^n\left(t\right)}$

is calculated using a numerical gradientwhere:

• • a=variable
  • c_t=compressibility (Pa⁻¹)
  • C=conductivity (m².m)
  • F_o,k=objective function (m³)
  • H=fracture height (m)
  • k_m=matrix permeability (m²)
  • L=maximum length of drainage (m)
  • m=real gas pseudo pressure (Pa/s)
  • m*=normalized real gas pseudo pressure (Pa/s)
  • m_f=real gas pseudo pressure in the fractures (Pa/s)
  • m_f*=normalized real gas pseudo pressure in the fractures (Pa/s)
  • m_m=real gas pseudo pressure in the matrix (Pa/s)
  • m_m*=normalized real gas pseudo pressure in the matrix (Pa/s)
  • m_—0=initial real gas pseudo pressure in the matrix (Pa/s)
  • M=molar mass (kg/mol)
  • P_m=matrix pressure (Pa)
  • P_m0=initial matrix pressure (Pa)
  • P_f=initial fracture pressure (Pa)
  • P_LB=low base pressure (Pa)
  • Q_tot=total flow rate from a matrix block into a fracture (m³/s)
  • Q_mf=local flow rate from a matrix block into a fracture (m²/s)
  • Q_f=flow rate inside the matrix (m³/s)
  • R=universal gas constant (J/mol/K)
  • t=time (s)
  • t_o,k=initial production time (s)
  • T=Temperature (K)
  • x_f=coordinate in the fracture (m)
  • x_m=coordinate in the matrix (m)
  • Z=volume factor
  • μ=viscosity (Pa.s)
  • φ_m=porosity
  • ρ=reservoir fluid density (kg/m³)
  • γ=variable

With t_o,k(t) known, the pressure distribution may then be calculated as follows:

$\rho Q_{\mathrm{mf}}\left(x_f,t\right)=\left(m_f-m_{\mathrm{m\_o}}\right)\frac{{\mathrm{MHk}}_m}{2\mathrm{RT}\sqrt{\pi a}}\sum _{k=1}^2\frac{1}{\sqrt{\left(t-t_{0,k}\right)}}\sum _{n=0}^{\infty }{\left(-1\right)}^n\left({}^{-{\left(\frac{\left(n+1\right)L_k}{\sqrt{a\left(t-t_{0,k}\right)}}\right)}^2}-{}^{-{\left(\frac{{\mathrm{nL}}_k}{\sqrt{a\left(t-t_{0,k}\right)}}\right)}^2}\right)$

• • This solution being linear in pressure, Eq. 2 may be integrated and solved.

$\begin{array}{cc}{m}_f\left(x_f\right)=m_{\mathrm{m\_o}}+\left(\frac{m_{f*,i}{}^{-\gamma L_f}-m_{f*,j}}{{}^{-\gamma L_f}-{}^{\gamma L_f}}\right){}^{\gamma x_f}+\left(\frac{m_{f*,j}-m_{f*,i}{}^{\gamma L_f}}{{}^{-\gamma L_f}-{}^{\gamma L_f}}\right){}^{-\gamma x_f}& \left(48\right)\\ {\gamma }^2\frac{k_m}{C\sqrt{\pi a}}\sum _{k=1}^2\frac{1}{\sqrt{\left(t-t_{0,k}\right)}}\sum _{n=0}^{\infty }{\left(-1\right)}^n\left({}^{-{\left(\frac{{\mathrm{nL}}_k}{\sqrt{a\left(t-t_{0,k}\right)}}\right)}^2}-{}^{-{\left(\frac{\left(n+1\right)L_k}{\sqrt{a\left(t-t_{0,k}\right)}}\right)}^2}\right)& \left(49\right)\end{array}$

Knowing the pressure distribution inside the network, the production rate can be calculated with Darcy's law.

The calculation may be implemented into various fracture networks without a constraint on the time step. In some cases, such as where production is done at a constant BHP and the conductivity is high, the flow from the matrix may be based on the assumption that the pressure inside the fracture stays constant. But in reality, only a fraction of the fracture branches of the network may have high conductivities. Examples of calculations are provided in Cipolla, C. L., Lolon, E. P., Mayerhofer, M. J., “Reservoir Modeling and Production Evaluation in shale-Gas Reservoirs”, SPE paper 13185 presented at the International Petroleum Technology Conference held in Doha, Qatar, Dec. 7, 2009.

Validity of the Analytical Solution

The analytical solution may be validated 427 by analyzing the solution to determine its application in a given formation. To study the validity of the analytical model for different values of fracture conductivity, the evolution of the pressure and production at a single fracture branch of a complex network may be analyzed. This study may consist of two sets of equally spaced parallel fractures as shown in FIG. 10. This figure describes a single branch 1070 in a DFN 1072 about a wellbore 1074 that will be analyzed. A matrix block 1048 of the DFN 1072 is depicted as having a flow rate 1076 from the matrix block 1048 to the fracture branch 1070.

For a high conductivity (finite) in the fracture network (e.g., about 2,500 mD.ft (762 mD.m)) in a reservoir of about 0.0001 mD, the BHP is almost instantaneously diffused inside the network and from there the variation of pressure inside the DFN can be neglected, compared to the pressure differential between the initial reservoir pressure and the BHP.

FIGS. 11.1 and 11.2 are three dimensional graphs 1100.1 and 1100.2 depicting reservoir pressure P (z-axis) versus distance x (m) (x-axis) and distance y (m) (y-axis) for 1 and 365 days, respectively. This figure depicts pressure of the DFN and an initial reservoir pressure 1178 at two different times of production for a high-conductivity DFN. These and other depictions may be provided. The production operation may be adjusted based on the production estimates.

As shown in FIG. 12, the pressure inside a selected fracture branch (such as the branch 1070 of FIG. 10) can be considered constant during ten years of production. This figure depicts a graph 1200 of pressure (P_m,o-P_f)(left y-axis) and initial time T in days (right y-axis) over time t in days (e.g., during three years of production) (x-axis) in a case of high conductivities. Generated lines for pressure 1280 and time delay 1281 are nearly flat.

The consequence of this almost constant pressure in the DFN is shown in FIG. 13 where the cumulated volume produced (from the matrix block into fracture branch (e.g., 1048 to 1070 in FIG. 10) converges towards the maximum volume recoverable as defined by mass balance (or initial volume of gas in place). FIG. 13 is a graph 1300 depicting cumulated production P (y-axis) versus time t (x-axis), resulting in a production curve 1384 that reaches toward a maximum recoverable volume 1382. This figure depicts cumulated production from a fracture branch versus time in case of high-conductivities. Because we are considering compressible fluids, in this example, the measure of volume may be done at surface conditions. This convergence indicates that the analytical solution verifies mass balance in this case of a high-conductivity DFN. The same analysis for a low-conductivity (or finite conductivity) DFN (50 mD.ft (15.24 mD.ft)) may have a different conclusion.

As shown in FIGS. 14.1 and 14.2, the pressure in the DFN may vary compared to the pressure range of the problem (e.g., BHP, initial reservoir pressure, etc.) This figure depicts pressure inside the DFN at two different times of production for low-conductivity DFN. FIGS. 14.1 and 14.2 are three dimensional graphs 1400.1 and 1400.2 depicting reservoir pressure P (z-axis) versus distance x (x-axis) and distance y (y-axis) for 1 and 365 days, respectively. This figure depicts pressure inside the DFN at two different times of production for a high-conductivity DFN. An initial reservoir pressure 1478 and a pressure of the DFN 1435 are also depicted.

This pressure variation may be seen on the pressure recorded in the fracture branch versus time as shown in FIG. 15. As shown in FIG. 15, where the pressure inside a selected fracture branch (such as the branch 1070 of FIG. 10) can be considered constant during ten years of production. This figure depicts a graph 1500 of pressure (P_m,o-P_f) (left y-axis) and time delay T in days (right y-axis) over time t in days (e.g., during three years of production) (x-axis) in a case of low conductivities (infinite). Generated lines for normalized pressure 1580 and time delay 1581 are nearly flat. Variation of boundary condition 1584 is also depicted.

This variation of the pressure in the DFN means that the assumption of constant pressure boundary condition in the analytical solution may require further analysis to confirm validity. The consequence is that the calculated flow rate from the matrix may be underestimated and the mass balance may be incorrect as illustrated in FIG. 16. FIG. 16 is a graph 1600 depicting cumulated production P (y-axis) versus time t (x-axis), resulting in a production curve 1684 that reaches toward a maximum recoverable volume 1682. This figure depicts cumulated production from a fracture branch versus time in case of low conductivities. An error 1686 between the production curve 1684 and the maximum recoverable volume 1682 is also depicted.

The low diffusivity in the fracture network may result in a “delay” in the production of the block depending on how far (or how connected) it is from the wellbore. This observation is a starting point for the method to extend the validity of the analytical solution to low conductivity fractures.

Extension of the Validity of the Analytical Solution

The validity of this analytical solution may be extended to modify the “initial” time t₀ (or t_o,k(t)) in such a way that the volume produced so far from the matrix block is equal to the volume that would have been produced according to the analytical solution under current pressure conditions inside the DFN. By doing this search at every time step and on each side of each fracture branch, the analytical solution is forced to satisfy mass balance. The search for the t₀ begins by defining the objective function F to minimize.

$\begin{array}{cc}{F}_k\left(t,t_{0,k}\left(t\right)\right)=M_{\mathrm{tot},k}\left(t\right)-{\int }_0^t{\int }_0^{L_f}\rho Q_{\mathrm{mf},k}\left(\tau ,P_f\left(x_f,t\right),t_{0,k}\left(t\right)\right)\tau x_f=0& \left(50\right)\end{array}$

M_tot is the volume produced at time t from the matrix block on the side k of the fracture branch. It is compared to the integration of the flow rate from the matrix over the length of the fracture branch and from the initial time t_0,k to t. The search for t_0,k such that F equals to zero, the iterative algorithm of Newton-Raphson as described in Eq. 51 may be used.

$\begin{array}{cc}{t}_{0,k}^{n+1}\left(t\right)=t_{0,k}^n\left(t\right)-\frac{F_k^n\left(t,t_{0,k}^n\left(t\right)\right)}{F_k^{\prime ,n}\left(t,t_{0,k}^n\left(t\right)\right)}\mathrm{with}F_k^{\prime ,n}\left(t,t_{0,k}^n\left(t\right)\right)=\frac{\partial F_k^n\left(t,t_{0,k}^n\left(t\right)\right)}{\partial t_{0,k}^n\left(t\right)}& \left(51\right)\end{array}$

The derivative of the function F_0,k is calculated by a numerical gradient. If t_0,k meets its time boundaries the optimization uses the bisection method. This optimization algorithm is very efficient because the solution from the previous time step is used as the initial guess for the next iteration. From a numerical point of view, the calculation of the approximated volume requires integration in time, which is the most CPU intensive part of the simulation. The optimization algorithm is applied for each side of each branch with minimal dependencies between the variables, making this part of the algorithm a candidate for parallel computing.

To illustrate the mechanism behind this approach, the analysis above may be used on a single fracture branch of a DFN with low conductivity (or finite conductivity) (50 mD.ft (15.24 mD.ft)). This pressure variation may be seen on the pressure recorded in the fracture branch versus time as shown in FIG. 17. As shown in FIG. 17, where the pressure inside a selected fracture branch (such as the branch 1070 of FIG. 10) can be considered constant during ten years of production. This figure depicts a graph 1700 of normalized pressure (P_m,o-P_f) (left y-axis) and time delay T in days (right y-axis) over time t in days (e.g., during three years of production) (x-axis) in a case of low conductivities. The resulting lines for normalized pressure 1780 and time delay 1781 incline.

FIG. 17 also shows the computed pressure inside the fracture and the initial time t_0,k updated with the proposed method. The increase of t_0,k with time may be necessary to sustain the flow rate from the matrix and the cumulated production as shown in FIG. 18. FIG. 18 is a graph 1800 depicting cumulated production (y-axis) versus time (x-axis), resulting in a production curve 1884 that reaches toward a maximum recoverable volume 1882. This figure depicts cumulated production from a fracture branch versus time in case of low conductivities.

This figure indicates that the method reduces the error on the mass balance, because the cumulated production converges close to the maximum recoverable volume in FIG. 18 in comparison to FIG. 16, thereby indicating that the validity of the method may extend with the analytical solution.

FIG. 19 is a chart 1900 depicting the distribution of pressure P and the initial time delay T as calculated by the algorithm on the entire fracture network at different time steps t₁ (1 day), t₂ (200 days) and t₃ (3 years). The chart includes DFNs 1935.1, 1935.2 and 1935.3 for pressure and DFNs 1935.4, 1935.5 and 1935.6 for time delay at the time steps t₁, t₂ and t₃, respectively. This figure shows pressure and initial time (or “delay’) in the reservoir at different times of production. The “pressure” column shows the pressure inside the reservoir blocks and the pressure inside the fracture network. The “initial time” T column shows the initial time for each block calculated by the algorithm.

The analysis above may be performed using an unconventional production model (UPM). To illustrate the capabilities of the UPM, simulations have been compared with those from a commercial reservoir simulator. The comparison is done with two different fracture geometries: a single bi-wing and a “wire-mesh” fracture network.

In an example involving a single bi-wing fracture, the hydraulic fracture is a single symmetrical fracture with a half-length of 1263 ft (384.96 m) and a fracture height of 98.4 ft (19.99 m). The permeability of the reservoir is 0.0001 mD with a porosity of 8%, the initial reservoir pressure is 4000 psi (281.29 kg/cm) and the bottom-hole pressure is 1000 psi (70.32 kg/cm). In this example, the volume factor Z and the gas viscosity were constant and equal to 1 and 0.02 cP, respectively.

FIG. 20 is a comparison of the simulated cumulated production between the reservoir simulation and the UPM, for different fracture conductivities varying between 0.005 and 5000 mD.ft (1524 mD.m), and for a bi-wing fracture. FIG. 20 shows that the greater the distance from the perforations (center of the grid), the smaller the initial time.

The FIG. 20 is a graph 2000 of cumulated production at surface condition (y-axis) versus time t (x-axis). This figure depicts a validation by comparison with a reservoir simulator. The resulting solid lines 2088.1-2088.7 and resulting dashed lines 2089.1-2089.7 show production based on reservoir simulator and the production model, respectively, at various locations. This graph 2000 indicates that the greater the distance is from the perforations, the longer it takes for the BHP to diffuse up to that location.

For a wire-mesh fracture network, this case represents a complex fracture network made up of 13 identical fractures in each orthogonal direction with a vertical well in the middle. In this example, the permeability of the reservoir is about 0.001 mD with a porosity of about 8%, the initial reservoir pressure is about 4000 psi (281.29 kg/cm) and the bottom-hole pressure is 1000 psi (70.32 kg/cm). As also shown in this example, the volume factor Z and the gas viscosity were constant and equal to 1 and 0.02 cP, respectively.

FIGS. 21.1 and 21.2 provide various visualizations of a DFN performed by various simulators. This figure show a reservoir and DFN used in a comparison between simulations done with a commercial reservoir simulator and the UPM. FIG. 21.1 shows an example of DFNs 2135.1 and 2135.2 as depicted by a reservoir simulator, such as ECLIPSE™. FIG. 21.2 shows a DFN 2135.3 generated using the UPM. Each of the DFNs depicted may be the same DFN resulting in the different images as shown.

FIGS. 22-24 compare results generated by a reservoir simulator and the UPM in examples where conductivities of DFN may vary. FIG. 22 is a comparison of the simulated cumulated production between the reservoir simulation and the UPM, for different fracture conductivities varying between 0.082 mD.ft (24.99 mD.mm) to about 8200 mD.ft (2499.36 mD.m), and for a bi-wing fracture.

FIG. 22 shows that the greater the distance from the perforations (center of the grid), the smaller the initial time. The FIG. 22 is a graph 2200 of cumulated production at surface condition (y-axis) versus time t (x-axis). This figure depicts a validation by comparison with a reservoir simulator. The resulting solid lines 2288.1-2288.6 and resulting dashed lines 2289.1-2289.6 show production based on reservoir simulator and the UPM, respectively, at various locations. This graph 2200 indicates that the greater the distance is from the perforations, the longer it takes for the BHP to diffuse up to that location.

As used herein, UPM without “delay” means that the UPM simulator uses the analytical part of the model, with a constant initial time equal to 0. When the fracture conductivity is increased, the difference between the reservoir simulator and the UPM simulation without “delay” may be reduced.

These comparisons show reasonably good agreement between the two simulators, in particular in the case of low conductivity where the algorithm to update the initial time plays a major role. To illustrate the importance of the algorithm correcting the initial time, FIGS. 23.1 and 23.2 compare the simulation results for the case of fracture conductivity equal to 82 mD.ft (24.99 mD.m).

FIG. 23.1 is a graph 2300.1 depicting flow rate at surface conditions. Production (y-axis) is plotted versus time t (x-axis). The resulting lines 2390.1-290.3 depict the simulation generated by a reservoir simulator, the UPM and the UPM without delay, respectively. FIG. 23.2 is a graph 2300.2 depicting current production at surface conditions. Cumulated production P (y-axis) is plotted versus time t (x-axis). The resulting lines 2390.4-290.6 depict the simulation generated by a reservoir simulator, the UPM and the UPM without delay, respectively. These figures depict a comparison of the rate (FIG. 23.1) and cumulated production (FIG. 23.2) between a commercial reservoir simulator, the UPM and the UPM without “delay”.

It should be noted that in the development of any such actual embodiment, numerous implementation—specific decisions must be made to achieve the developer's specific goals, such as compliance with system related and business related constraints, which will vary from one implementation to another. Moreover, it will be appreciated that such a development effort might be complex and time consuming but would nevertheless be a routine undertaking for those of ordinary skill in the art having the benefit of this disclosure. In addition, the composition used/disclosed herein can also comprise some components other than those cited. In the summary and the detailed description, each numerical value should be read once as modified by the term “about” (unless already expressly so modified), and then read again as not so modified unless otherwise indicated in context. Also, in the summary and the detailed description, it should be understood that a concentration range listed or described as being useful, suitable, or the like, is intended that any and every concentration within the range, including the end points, is to be considered as having been stated. For example, “a range of from 1 to 10” is to be read as indicating each and every possible number along the continuum between about 1 and about 10. Thus, even if specific data points within the range, or even no data points within the range, are explicitly identified or refer to only a few specific, it is to be understood that inventors appreciate and understand that any and all data points within the range are to be considered to have been specified, and that inventors possessed knowledge of the entire range and all points within the range.

The statements made herein merely provide information related to the present disclosure and may not constitute prior art, and may describe some embodiments illustrating the disclosed subject matter. All references cited herein are incorporated by reference into the current application in the entirety.

The preceding description has been presented with reference to some embodiments. Persons skilled in the art and technology to which this disclosure pertains will appreciate that alterations and changes in the described structures and methods of operation can be practiced without meaningfully departing from the principle, and scope of this application. Accordingly, the foregoing description should not be read as pertaining only to the precise structures described and shown in the accompanying drawings, but rather should be read as consistent with and as support for the following claims, which are to have their fullest and fairest scope.

Although only a few example embodiments have been described in detail above, those skilled in the art will readily appreciate that many modifications are possible in the example embodiments without materially departing from the system and method for performing wellbore stimulation operations. Accordingly, all such modifications are intended to be included within the scope of this disclosure as defined in the following claims. In the claims, means-plus-function clauses are intended to cover the structures described herein as performing the recited function and not only structural equivalents, but also equivalent structures. Thus, although a nail and a screw may not be structural equivalents in that a nail employs a cylindrical surface to secure wooden parts together, whereas a screw employs a helical surface, in the environment of fastening wooden parts, a nail and a screw may be equivalent structures. It is the express intention of the applicant not to invoke 35 U.S.C. §112, paragraph 6 for any limitations of any of the claims herein, except for those in which the claim expressly uses the words ‘means for’ together with an associated function.

Claims

1. A method of performing a production operation about a wellbore penetrating a subterranean formation, the subterranean formation having a plurality of fractures thereabout, the method comprising: generating a flow rate through a discrete fracture network, the discrete fracture network extrapolated from a hydraulic fracture network defined by the plurality of fractures in the subterranean formation, the discrete fracture network comprising a plurality of fracture branches with intersections therebetween and a plurality of matrix blocks;generating a pressure profile of the discrete fracture network for an initial time based on the flow rate; andgenerating a production rate based on the pressure profile.

2. The method of claim 1, wherein the generating the flow rate comprises generating the flow rate from one of the plurality of matrix blocks to one of the plurality of fracture branches.

3. The method of claim 1, wherein the generating the flow rate comprises generating the flow rate inside at least one of the plurality of fractures.

4. The method of claim 1, wherein the generating the flow rate comprises generating the flow rate from one of the plurality of matrix blocks to one of the fracture branches.

5. The method of claim 1, wherein the generating flow rate comprises generating flow rate inside the one of the fracture branches between two of the intersections of the discrete fracture network.

6. The method of claim 1, wherein the generating the flow rate comprises generating the flow rate inside the one of the plurality of fracture branches at the intersections of the discrete fracture network.

7. The method of claim 1, further comprising determining mass balance at the intersection between two of the plurality of fracture branches.

8. The method of claim 1, wherein the generating the pressure profile comprises generating the pressure profile using Darcy's law.

9. The method of claim 1, wherein the generating the pressure profile is unconstrained by a time step.

10. The method of claim 1, further comprising defining a time function of fluid flow through the matrix block, the time function having the initial time.

11. The method of claim 10, further comprising updating the time function of fluid flow through the matrix block.

12. The method of claim 1, further comprising updating the production rate at a plurality of time steps.

13. The method of claim 1, further comprising updating the production rate for the plurality of fracture branches.

14. The method of claim 1, further comprising accounting for delays in production of each of the plurality of matrix blocks by updating the initial time such that an actual mass produced from each of the plurality of matrix blocks equals the mass if the current pressure conditions in an adjacent one of the plurality of fracture branches would have been constant in time from the updated initial time.

15. The method of claim 1, further comprising validating the production rates.

16. The method of claim 15, wherein the validating comprises comparing the production rates with production rates generated by a reservoir simulator.

17. The method of claim 15, wherein the validating is performed for the discrete fracture network having high conductivity, low conductivity, bi-wing fractures, wire-mesh fractures, time delay, and combinations thereof.

18. The method of claim 15, wherein the validating comprises modifying the initial time so that a volume produced from the plurality of matrix blocks over time for each of the plurality of fracture branches satisfies mass balance.

19. A method of performing an oilfield operation about a wellbore penetrating a subterranean formation, the subterranean formation having a reservoir therein, the method comprising: performing a fracture operation, the fracture operation comprising generating fractures about the wellbore, the fractures defining a hydraulic fracture network about the wellbore;generating a discrete fracture network about the wellbore by extrapolating fracture data from the hydraulic fracture network, the discrete fracture network comprising a plurality of fracture branches with intersections therebetween and a plurality of matrix blocks;generating a depth of drainage through the discrete fracture network;defining at least one production parameter; andperforming a production operation to produce fluids from the subterranean formation based on the depth of drainage and the at least one production parameter.

20. The method of claim 19, further comprising measuring downhole data about the wellbore.

21. The method of claim 19, wherein the performing the fracture operation comprises stimulating production from the wellbore by injecting fluid into the subterranean formation.

22. The method of claim 19, wherein the performing the fracture operation comprises simulating the performing the fracture operation.

23. The method of claim 19, wherein the discrete fracture network considers an averaged value for at least one fracture property at each of the plurality of fracture branches.

24. The method of claim 23, wherein the at least one fracture property comprises spatial coordinates at a fracture branch extremity, conductivity, averaged conductivity, height, averaged height, reservoir pressure, averaged reservoir pressure at a fracture branch location, permeability, averaged reservoir permeability at the fracture branch location and combinations thereof.

25. The method of claim 19, wherein the generating the depth of drainage comprises evaluating the depth of drainage through the plurality of matrix blocks of the discrete fracture network.

26. The method of claim 19, wherein the generating the depth of drainage comprises generating the depth of drainage for each of the plurality of matrix blocks based on an approximation of linear flow through the plurality of matrix blocks.

27. The method of claim 19, wherein the generating the depth of drainage comprises automatically evaluating the depth of drainage of the plurality of matrix blocks to be depleted in front of each of the plurality of fracture branches and accounting for a volume to deplete for each of the plurality of matrix blocks.

28. The method of claim 19, wherein the at least one production parameter comprises bottom hole pressure, reservoir fluid viscosity at reservoir conditions, reservoir fluid compressibility at reservoir conditions, duration over which production is to be simulated, and combinations thereof.

29. The method of claim 19, wherein the performing the production operation comprises positioning tubing in the wellbore and transporting fluids from the reservoir to a surface location.

30. The method of claim 19, wherein the performing the production operation comprises estimating a production rate from the wellbore by simulating the production of fluid from the wellbore.

31. The method of claim 30, wherein the performing a production operation comprises visualizing the production rate.

32. The method of claim 30, further comprising adjusting the performing based on the estimated production rate.

33. The method of claim 19, wherein the performing the production operation is based on a range of fracture parameters.

34. A method of performing an oilfield operation about a wellbore penetrating a subterranean formation, the subterranean formation having a reservoir therein, the method comprising: stimulating the wellbore by injecting fluid into the subterranean formation such that fractures are generated about the wellbore;measuring the fractures and defining a hydraulic fracture network based on the measured fractures;generating a discrete fracture network about the wellbore by extrapolating fracture data from the hydraulic fracture network, the discrete fracture network comprising a plurality of fracture branches with intersections therebetween and a plurality of matrix blocks;generating a depth of drainage through the discrete fracture network;defining at least one production parameter; andestimating a production rate over time based on the depth of drainage and the at least one production parameter; andproducing fluids from the subterranean formation based on the estimated production rate.