2. TASK PROGRESS

2.1 Active Tasks
2.2 Task Progress
2.2.1 Task 1.1.2: Data Updates
2.2.2 Task 1.2.1: 3D Hierarchical Fracture Model
2.2.3 Task 1.2.2 Hierarchical Fracture Model Verification
2.2.4 Task 3.1.1: Linkage to Reservoir Models (Software Development)
2.2.5 Task 3.2.1: MS Windows 95 Analysis System
2.2.6 Task 3.2.3: Software Linking
2.2.7 Task 4.1.2: Well Testing Data Acquisition
2.2.8 Task 5.1.2: WWW Site Updates
2.2.9 Task 5.2.1 Progress Reports
2.2.10 Task 5.2.2: Research Reports
2.2.11 Task 5.2.3: Presentations
2.2.12 Task 6: Management

2. TASK PROGRESS

2.1 Active Tasks

The following tasks were active during the quarter:

• Task 1.1.2 Data Updates
• Task 1.2.1 3D Hierarchical Fracture Model
• Task 1.2.2 3D Hierarchical Model Verification
• Task 3.1.1 Linkage to Reservoir Models/StrataModel
• Task 3.2.1 MS Windows 95 Analysis System
• Task 3.2.3 Software Linking/ECLipse
• Task 4.1.2 Well Testing Data Acquisition
• Task 5.1.2 WWW Site Updates
• Task 5.2.1 Progress Reports
• Task 5.2.3 Presentations
• Task 6 Management

2.2 Task Progress

This section describes progress during the quarter for each of the active tasks.

2.2.1 Task 1.1.2: Data Updates

During the quarter, Marathon Oil Company collected production and field testing data from the project study site and provided this data to Golder Associates.

Well performance histories were provided as monthly oil, gas, and water rates from each well within the two window areas for each of over 120 wells in the project study site. The production was allocated based on monthly individual well test data, actual producing time during the month for each well, and actual metered production from groups of wells throughout the month.

The individual wells and group total for the Tract 49 area reflect the performance of wells responding to the TAGS process. Currently these wells are recovering the water condensate from the injected steam but not at a significantly altered water chemistry than the previous production water (the water is equilibrating with the carbonate formation and acid gases of the reservoir (CO₂ and H₂S) prior to being produced).

Well test data is provided (on all active producers within the window areas) as actual metered test periods of three phase production (typically 24 hr. test of oil, gas, and water rates). This data provides the most accurate indication of water-oil ratio and gas-oil ratio (WOR and GOR) changes that result from the oil column moving down to allow deeper discrete features to flow oil into the wellbores, or to see the loss of oil rate as the local oil column settles into an elevation where specific wells have no good discrete connection to the oil-filled fracture network.

For comparison with, and assessment of, the production well test data, static observation well data is provided for several dozen locations in or near the window areas. This critical data provides the location and rate of movement of the gas-oil and water-oil contacts (GOC and WOC) within the fracture network. Oil is only able to flow through network connections within the boundaries of the GOC and WOC. Network connection beyond either fluid contact are not oil phase connections.

Pressure measurements are provided in these observation wells to monitor changes in the energy level and to confirm the multi-phase "piezometric surface" or isobaric surface at datums within each phase (gas, oil, and water). As the oil column drops at 14 feet per year, the discrete flow connections continuously change within the oil saturated fracture network and to individual wellbores withdrawing the oil. As connected features empty of oil, well performance deteriorates until the oil moves into a deeper feature connected to the wellbore.

A database of fluid samples and flowing temperatures is being collected, analysed, and organized into a database for analysis of TAGS response. This database is expected to be added to the project Web site during the next quarter. Review of this data indicates a background (baseline) scatter of compositions in each phase with a few anomalous points.

Marathon has also provided multi-well pressure interference and chemical tracer test results for the Tract 17 area for evaluation with the continuing performance, pressure, and contact data.

The data provided by Marathon during the quarter was assembled and is being incorporated into the project World-Wide-Web (WWW) site (Task 5.1.2). Figures 2-1, 2-2, 2-3, and 2-4 present examples of the data provided by Marathon during the quarter.

2.2.2 Task 1.2.1: 3D Hierarchical Fracture Model

During the quarter, MIT continued the research in two major areas: 1) further theoretical development of the 3D hierarchical fracture model (HFM) (Ivanova 1995) and 2) numerical simulation of the Yates reservoir fracture network. The theoretical development included modeling of fracture intensity, and programming of algorithms applicable to fracture networks related to folds. The numerical simulations for the Yates field focused on the fracture network associated with the asymmetric dome structure.

Theoretical development during the quarter focused on

• stress field orientation,
• stress field variation, and
• intensity and fracture size.

2.2.2.1 Stress Field Orientation

In the HFM, orientations of potential fracture planes are represented by the primary stochastic process: a homogeneous, anisotropic, Poisson plane network. The mean orientation of a fracture set is specified in polar coordinates in terms of an azimuth Q and a latitude F (Figure 2-5a) in a global frame of reference (OXYZ). In the frame of reference of a fracture set (Oxyz in Figure 2-5a) is a plane is defined by the equation:

(Equation 2-1)

where d is the orthogonal distance from the origin to the plane, and q and f are the latitude and azimuth of the normal vector (pole Oz'), respectively. A homogeneous Poisson plane network of intensity m corresponds to a Poisson point process in:

(Equation 2-2)

with non-homogeneous intensity function of the type:

(Equation 2-3)

where f_q,f is the joint PDF of q and f and m is a positive constant. The ordered distances from an arbitrary point to the planes of a Poisson network with intensity m define a Poisson point process on a line with intensity 2m (Miles 1969; Veneziano 1978). This property is used in the algorithm for generation of random planes described by equation (1): directional pairs (q, f) fit a spherical PDF (Dershowitz 1979), and distances d from the origin of the global frame of reference follow a Poisson point process of intensity 2m.

The plane network is generated in a modeling volume that represents the fractured rock. Figure 2-5b illustrates the modeling volume, defined by four vertical planes at X=X_m, X=- X_m, Y= Y_m, Y=- Y_m, one horizontal plane at Z=0, and a curved top surface, described by a quadratic function:

(Equation 2-4)

The top surface can represent any geologic structure that limits the extent of propagation of the fracture network, for example, a topographic surface, a fold, or a bedding plane.

2.2.2.2 Stress Field Variation

The rotation/translation procedure of the tertiary stochastic process represents possible deviations of fracture orientations from the regional directions due to local variations of the stress field. In the tertiary stochastic process, translation is performed in the frame of reference of the fracture plane (O'x'y'z' in Figure 2-5a). Translation is accomplished by assigning a non-zero coordinate z'i to the center and the vertices of a polygon i (hence to the entire polygon). The algorithm, currently implemented in the model, allows for translation of a polygon at a maximum distance:

(Equation 2-5)

where Re is the radius of a circle with equivalent area, and E[R_e] is the expected value of that radius. Thus larger polygons are shifted closer to their original position than smaller polygons. C in Equation 2-5 is a coefficient of fracture co-planarity: small C indicates more co-planarity of fractures in the same set. Figure 2-6 illustrates the translation process in 2-D. The fracture systems shown on horizontal outcrops in (b) and (c) are statistically the same (in terms of P₃₂, mean size, orientation, etc.); only the translation coefficient is different. The figure illustrates the effect of increasing C to decrease the co-planarity of the fractures.

The random rotation of polygons in the tertiary process is used for cases when the conditions modify the stress field so much that locally the most likely fracture orientation deviates from the most likely orientation defined by the general stress field. The algorithm, currently implemented in the model, checks the orientations of synthetic polygons against a specified relation of the real fractures to the orientation of geologic structures, bedding planes, or other features that may influence the direction of fracture propagation. Figure 2-7 illustrates the rotation process. In (a) the system is generated with intensity P₃₂ and a PDF related to the orientation of the general stress field. In (b), after rotation, fracture strikes are concentric or radial to a circular dome with apex in the center of the rectangular area. The rotation process preserves the fracture intensity P₃₂ of the fracture set.

2.2.2.3 Fracture Intensity and Fracture Size

In the HFM, fracture intensity is defined as the cumulative fractured surface area per unit volume of rock:

(Equation 2-6)

where Af,i is the area of a fracture inside the volume V. P₃₂ is a good measure of fracture intensity because (1) it represents the two-dimensional nature of fractures, and (2) it does not depend on the size of the sampled region as long as it is representative for the fracture network (Dershowitz and Herda, 1992). In the 3D hierarchical model, fracture intensity is modeled by the combined primary and secondary stochastic processes. Some characteristics of the Poisson point, line and plane processes, previously published (Miles 1969; Miles 1973; Veneziano 1978), or established herewith via numerical simulations, allow for the development of a straightforward procedure for generation of populations of polygons-fractures with desired intensity of the type defined in Equation 2-6.

A random subdivision of the planes generated by the primary process into a fractured region and its complementary region of intact rock constitutes the basis of the secondary process. This subdivision is accomplished by a Poisson line tessellation on every plane and a non-homogeneous process of marking the so-created polygons as fractured or intact rock.

If A is the polygonal intersection between the fracture plane and the modeling volume V and O'x'y' is the local 2D frame of reference on the plane (Figure 2-5a), a line from the Poisson network is defined as:

(Equation 2-7)

where a is an angle measured counterclockwise from the axis O'x', and D is the orthogonal distance from the origin O' to the line.

A homogeneous Poisson line network with intensity l corresponds to a Poisson point process in the region:

(Equation 2-8)

with intensity function of the points of the type:

(Equation 2-9)

where l is a positive constant and f_a(a) is the PDF of a in [0,2p] (uniform in the model). The ordered distances from an arbitrary point on the plane to the lines of a Poisson line network with intensity l form a Poisson point process with intensity 2l. The points of intersection of all lines determine the vertices of polygonal shapes.

Table 2-1 summarizes first and second moments of the polygons created by a line tessellation of intensity l.

Table 2-1 Statistics of polygons by a Poisson line tessellation with
intensity l (after Miles 1973, and Veneziano 1978)

The standard deviation of polygon areas can be calculated as:

(Equation 2-10)

In the model, the process of dividing each plane, A, into a fractured region, A1, and its complementary region of intact rock, A0, is non-homogeneous. The probability Pf of marking a polygon as fractured is calculated individually for every polygon:

P_f = P_f (size, shape, location, ...) (Equation 2-11)

A polygon, produced by the line tessellation, is retained as a possible fracture only if it has an appropriate shape. The mean and variance of the size of polygons that are potential fractures are different from that of the entire population (Table 2-1). Approximate First Order Second Moment analysis calculates the conditional mean and variance for polygons with a certain number of vertices as:

(Equation 2-12)

(Equation 2-13)

For example, the mean area of the polygons with five, six, and seven vertices is E[A/N=5]=2E[A], E[A/N=6]=3E[A], and E[A/N=7]=4E[A], respectively, where E[A]=p/l2.

The desired mean size of the fracture set members can be related to the intensity of the Poisson line network necessary to produce the polygons. Table 2-1 summarizes the results of numerical generations, performed in order to establish two correction coefficients:

• the ratio CA of the mean area E[Am] of "good polygons" to the mean area E[A] of all polygons produced by a tessellation with intensity l;
• the percentage g of cumulative area of "good polygons" per unit area.

Table 2-2 Calculated statistics of marked polygons created by a Poisson line tessellation

In the 3D hierarchical model, "good polygons", i.e., potential fractures, are marked according to three rules that reflect the shapes typical for natural fractures:

• the polygons have at least four vertices;
• all angles are at least p/3 (60 degrees);
• the polygons are not extensively elongated.

The results for every one of the five cases in Table 2-2 are based on at least 100 simulations of Poisson line tessellation on a plane with intensity calculated as:

(Equation 2-14)

where E[A] is the desired mean area of the total population of fractures. The population of polygons with appropriate shapes for the purpose of fracture network modeling has the following properties:

Mean area:

(Equation 2-15)

Standard deviation:

(Equation 2-16)

Cumulative area of "good polygons" per unit total area:

(Equation 2-17)

In order to obtain the expected value of fracture intensity P₃₂ (Equation 2-6), one needs to know not only g but also the expected value of the cumulative total area A_T,f that will be cut from a volume V by a Poisson network of planes with intensity m. For example, a plane at a random distance d cuts from a sphere with radius R a circle with expected area E[A_T]:

(Equation 2-18)

where pdf(d)=1/R is the probability density function of d. The expected cumulative area of potential fractures per unit volume is therefore:

(Equation 2-19)

where 2mR is the expected number of planes from the Poisson network that intersect the sphere. The expected cumulative area of polygons-potential fractures per unit volume does not depend on R, i.e., it is independent on the shape and size of the total volume in which the Poisson plane process is generated. It only depends on the intensity of the Poisson plane network m and on the rule by which polygons are marked as potential fractures. According to the rule which has been chosen here, g=0.4.

There are two ways to obtain a smaller percentage g* of the potentially fractured planes (i.e., a lower percent of coplanar fractures). One can either use the tertiary stochastic process, or mark the "good polygons" with a probability P such that:

(Equation 2-20)

where 0<P<1.

Then a Poisson plane network with intensity m*=m/P preserves the fracture intensity:

(Equation 2-21)

In the 3D model, polygons are marked with a different zone probability as a function of their distance to some specified cluster centers, for example, other fractures or local faults. A zone is defined by the maximum and minimum value of the distance from the face of a polygon-fracture to the face of another polygonal feature. For example, Figure 2-8 illustrates a fracture network with higher intensity in the vicinity of a fault. Figure 2-9 illustrates a HFM containing two fracture sets: one (a) defined in a layer enclosed between two bedding surfaces but not in the surrounding rock, and (b) defined with high intensity in the rock below a bedding plane than above it.

Within every zone j the zone probability P_j is constant, but that constant is generally different in different zones. Polygons in zone j are retained with probability P_j and discarded with probability 1- P_j. To obtain different fracture intensity within the same fracture set in different portions V₁, V₂, etc. of the volume V, one has to calculate the different probabilities P₁, P₂, etc. for polygon marking in the different regions. Volume V_i does not have to be continuous; it can be defined as a function of any continuous rock property. For example, V_i may be defined as "the regions where the fractures of the set are at a distance not more than L" from the fractures of a primary set or from a fault face. Alternatively, V_i can be defined as "the regions where the rock is dolomite" or "the regions where the porosity of the rock is not more than n%", etc.

2.2.3 Task 1.2.2 Hierarchical Fracture Model Verification

During the quarter, verification of the HFM was initiated through simulation of the well-scale fracture network in Tract 49. The well scale fracture model is developed in the context of larger and smaller scale fractures in Dershowitz et al. (1996, 1997). The model is developed within the context of a larger scale asymmetric domal structure.

The present day structure of the M-horizon in the Seven Rivers formation (most likely deposited as a horizontal layer in the Late Permian) has been used for inference of the fold geometry to which the fracture network of interest is related. Figure 2-10 illustrates the M-horizon and its approximation by a cubic surface of the type:

(Equation 2-22)

determined by a polynomial fit to the elevations of the M-horizon observed in logged wells throughout Tract 49.

Table 2-3 summarizes the derivation of parameters for the HFM simulation of Tract 49.

Table 2-3 Parameters of the 3D hierarchical model used in simulations
of the fracture network in Tract 49

Every numerical generation is performed in a fracture generation region of the type illustrated in Figure 2-5b where:

• The origin of the frame of reference is centered at one of the five logged wells in Tract 49 from which fracture FMI/FMS interpretation is available: YU3728, YU4007, YU4903, YU4959, or YU5127.
• Datum Z=0 is at elevation 700 ft above sea level.
• The upper bounding surface is assumed at about 100 ft above the San Andres formation top which is represented by a quadratic surface (Equation 2-4) fitted to elevations observed in logged wells.
• X_m = 400 ft and Y_m = 400 ft; thus the horizontal extent (800 ft) is approximately equal to the average vertical thickness between the datum and the top boundary.

Figures 2-11, 2-12, 2-13, and 2-14 illustrate the numerical simulations for well YU4007. In Figure 2-11 contours of the fold surface in the vicinity of the well are shown in (a), and fracture traces striking predominantly sub-parallel to the strike of the fold are shown in (b) on a horizontal cross-section. Figure 2-12 depicts a 3D view of the numerically generated fractures intersected by a hypothetical well. Figure 2-13 compares fracture strike diagrams of the borehole intersections of the synthetic network in three independent simulations to the strike rosette diagram of fractures observed in the available FMI/FMS log at well YU4007. The match of strike orientation distribution is very good. The difference in fracture number and spacing between the real and the synthetic fracture networks is due to the fact that fracture intensity observed in cores (lower than that observed in logs) has been used in the numerical simulations. Figure 2-14 compares dip distribution of the fracture network from log analysis at well YU4007 to the dip distribution from an HFM Simulation. The peak of the dip distribution of the numerical simulation (80-85 degrees) is shifted from that of the fractures observed in logs (about 65-70 degrees). This is again due to the assumption in the inference procedure that core observations are more reliable than log interpretations: fractures in cores have dips predominantly 80 degrees and higher.

Figure 2-15 compares strikes of HFM Model Fracture to FMI/FMS measurements from logs at wells YU3728, YU4903, YU4959, and YU5127 in Tract 49. In Figures 2-11, 2-12, 2-13, 2-14, and 2-15, only those numerically generated fractures that are intersected by simulated outcrop planes or boreholes are shown. The total number of numerically generated fractures in every simulation for a given well exceeds 100,000 due to the relatively high fracture intensity P₃₂ inferred from field data.

2.2.3.1 Outlook

The HFM model verification task will continue during the upcoming quarter with the following activities:

• Direct linking of StrataModel dolomite porosity data (Dershowitz et al. 1996) for fracture generation.
• Numerical generation of the fracture network for the Tract 17 study site. Tract 17 exhibits a more complex relationship to the field structure exists, possibly combined with fracture termination at shale layers and influence of field-scale faults/drape folds on fracture intensity.

The preliminary results indicate that the intermediate-scale fracture system in the Yates Field is likely to have a good connectivity in relation to the field dome. In Tract 49 the connectivity is mainly concentric to the dome (Figure 2-16). The 3D fracture network illustrated in Figure 2-16 shows only a portion of the fracture network, with intensity much lower than P₃₂=0.4 inferred from field data for Tract 49. It is important to note that the intensity of intermediate-scale fractures that may provide significant conduits for flow and transport in the Yates Field is probably lower than that of the entire network. In the inference of conductive fracture intensity from field data it is essential to consider the 3D distribution of fracture filling secondary calcite and the results from field water injection tests.

2.2.4 Task 3.1.1: Linkage to Reservoir Models (Software Development)

Prototype software was developed during the quarter to provide a linkage between FracMan discrete feature network. The prototype "StrataFrac" has two functions:

1. Convert the StrataModel *.RW file format into a file format to facilitate direct generation of fractures based on StrataModel grid cell values
2. Generate StrataModel *.RW grid cell values from 3D DFN geometries.

The flow chart for StrataFrac illustrated in Figure 2-17. The algorithm is described in Dershowitz et al. (1997). The format of the StrataModel *.RW file is illustrated in Table 2-4.

Table 2-4 StrataModel *.RW File Format Example

The algorithms implemented are described in Dershowitz et al. (1997). The StrataFrac user interface is illustrated in Figure 2-18.

2.2.4.1 StrataFrac User Interface

To execute StrataFrac, select the StrataFrac icon, and double click with the mouse. The general procedure for analysis is summarized in Table 2-5. The operation of the individual StrataFrac menu items are described in the next section. Navigation through StrataFrac is done using Microsoft Windows mouse conventions. In general, the left hand mouse button is used for making selections.

Table 2-5 StrataFrac Analysis Sequence: DFN_StrataModel

2.2.4.2 Command Summary

2.2.4.2.1 File Menu

New: Begin a new analysis, without closing the current analysis. Opens a window for the new analysis.

Open: Open fracture geometry (.FAB) and StrataModel (*.RW) files. The standard Microsoft Windows File Open menu syntax and assumptions are used.

Close: Close the current active analysis and all windows related to that analysis.

Save As: Saves fracture intensity field (.DAB) and StrataModel (*.RW) files

Exit: Exit StrataFrac. Warns the user if statistical reports have not been saved.

2.2.4.2.2 Edit Menu

Undo: Not currently implemented

Cut: Not currently implemented

Copy: Not currently implemented

Paste: Not currently implemented

Options: Define the transformation to be performed:

• type of transformation (DFN_StrataModel or StrataModel_DFN)

  transformation and rotation between DFN and StrataModel coordinate systems

  slab region of DFN or StrataModel to be used as the basis of the transformation

  For StrataModel_DFN transformation, the StrataModel variable name to be used to define the DFN "probability of fracture center" for a cell

  For DFN_StrataModel, the names of the StrataModel variables to contain the Is_Frac, N_Frac, and Intensity P₃₂, computed from the DFN

2.2.4.2.3 View Menu

Toolbar: Display or hide the toolbar display

Statusbar: Display or hide the line providing status information.

Fractures: Display or hide fractures.

StrataModel: Display or hide StrataModel cells

2.2.4.2.4 Window Menu

New Window: Open a new window for the current analysis

Cascade: Cascade currently open windows

Tile: Tile currently open windows

Arrange Icons: Neatly arrange icons for minimized windows

Close: Close the current window

Close All: Close all the windows in the current analysis.

2.2.4.2.5 Help Menu

Help Index: Not currently implemented

Using Help: Not currently implemented

About StrataFrac: StrataFrac QA, copyright, and license information.

2.2.4.3 StrataFrac Walk-through

This section provides a walk-through of the StrataFrac application (Figure 2-18). After launching StrataFrac, the user is confronted with the information and licensing screen. Upon accepting the license agreement by clicking on "OK", the user comes to the main window. Under the Files menu, the user should use the mouse to highlight the discrete fracture network (*.FAB) file and the StrataModel file (*.RW) to be processed.

Returning to the main menu, the user next selects Edit/Options, to define the way in which the interface is to be used to convert between DFN and StrataModel. Once these have been set, the user can continue to File/Save As to save the converted file.

2.2.5 Task 3.2.1: MS Windows 95 Analysis System

Flare and Fractal are still under development. The other components of the MS Windows 95 analysis system are complete, operational, and are now available on the World Wide Web. All five programs will undergo further improvements and should be considered Beta 1.0 versions.

The status of active development of software is summarized in Table 2-6.

Table 2-6 MS Windows 95 Analysis System

2.2.6 Task 3.2.3: Software Linking

The linkage between DFN models and the ECLipse reservoir simulator was further developed during the quarter.

Reservoir simulation for fractured reservoirs can be significantly more challenging than it is for conventional porous media reservoirs. Issues of connectivity and heterogeneity which play a key role in fractured reservoirs can frequently be ignored in simplified equivalent porous medium (EPM) approaches for conventional reservoirs. The most popular approach for analysis of fractured reservoirs is to extend the EPM approach of porous reservoirs to fractured reservoirs by adding a second interacting continuum. This dual porosity continuum (DP) approach addresses only the dual porosity nature of fractured reservoirs, but completely ignores the connectivity and heterogeneity issues.

However, despite the significant simplifications made in DP models regarding the geometry of the fracture network, they are still the leading approach in practice. As a result, there is a need for development of linkages between DP and DFN approaches. These linkages will allow the DP approach to take advantage of some of the features of the DFN approach, without requiring reservoir engineers to be extensively retrained. It also allows the analyst to maintain many of the advantages of DP simulators such as ECLipse, including:

• computational efficiency at reservoir scales
• multiphase flow modeling
• coupled process modeling
• highly developed EOR simulation

This task develops an approach for integration of the realism of complex fracture system geometries, as captured by DFN models, with the complexity of flow simulations offered by DP continuum models. The approach adopted is to derive the continuum parameters with the help of DFN models. This approach requires that DFN models are used to provide the parameters needed for reservoir simulation, on a cell by cell basis. In We have developed approaches to derive the following DP model parameters from DFN models:

• optimum grid cell size
• fracture system porosity
• directional fracture system permeabilities
• matrix-fracture interaction factors (spacing, _-Factor)
• permeability barriers
• marking of active and inactive cells
• well permeability thickness Kh
• well productivity/injectivity index

The analysis approaches for each of these parameters are described below in Section 2.2.6.1. A software package, "FracLips" is under development to provide these features.

It is important to note that while the FracLips approaches can generally be applied to transform DFN models to ECLipse based DP models, this does not necessarily guarantee that the resulting ECLipse DP model will encapsulate the true connectivity and heterogeneity of the fractured reservoir. This can frequently be achieved only by applying the DFN approach directly. Further, it only makes sense to transform DFN data to DP models, if the rock mass can be treated as an EPM continuum. This is mainly a function of scale. The cell size analysis shows whether this is the case or not. The DFN approach does NOT need any upscaling or downscaling because it uses the actual grid cell size to compute permeabilities.

The transformation of DFN models to continuum parameters can be carried out as:

• DFN_Homogeneous Continuum, in which the DFN is used to derive a "representative" grid cell property to be applied to every grid cell
• DFN_Heterogeneous Continuum, in which the DFN is used to derive the properties of specific grid cells, where each grid cell can have different values
• DFN_Stochastic, Homogeneous Continuum, in which the DFN is used to derive distributions of the properties of a "representative" grid cell.
• DFN_Stochastic, Heterogeneous Continuum, in which multiple realizations of the DFN are used to derive multiple realizations for the properties of specific grid cells, where each grid cell can have different values.

Which of those options is used depends on how the continuum simulations are to be carried out.

The DFN approach is based on the stochastic modeling concept. Hence every realization of the discrete fracture network will produce different results. Multiple realizations need to be carried out and the results of the individual realizations have to be post-processed to get final results. Currently most continuum models are as single `best guess' realizations, i.e., only one geological model is set up which is then history matched. In this case the model represents the best guess for the fracture system parameters based on multiple DFN model realizations (e.g., mean porosity from all DFN model realizations).

In the future, more reservoir studies will use probabilistic simulation methods (Eiben et al., 1996; Lia et al., 1997). In this case multiple realizations of the continuum model are used to honor uncertainties in the input parameters. For those studies FracLips will be able to sample the parameters of a DFN realization onto an existing ECLipse simulation grid. The continuum model parameters will only be valid for the unique fracture network in this realization. The sampling and subsequent modeling efforts be repeated for multiple realizations in order to provide a valid description of the reservoir.

2.2.6.1 Approach for Calculation of DP Continuum Parameters

2.2.6.1.1 Optimum Grid Size

From a discrete fracture network perspective, the optimum grid size would be one at which each grid cell is hydraulically connected to each of its surrounding grid cells. This is illustrated in Figure 2-19. At scales smaller than the optimal grid size, a significant proportion of the grid cells are not connected to their surrounding grid cells. The DP approach assumes connectivity between each grid cell and its six immediate neighbors. As a result, unless the scale is large enough to provide this level of connectivity, the hydraulic behavior of the DFN and DP representations are quite different.

To ensure that the grid cell size in the DP approach is large enough, FracLips sets up a series of grid cells systems at a range of sizes specified by the user. FracLips then uses graph theory searches to determine the number of connected grid cells per grid cell for each discretization. An example of this analysis is presented in Figure 2-20.

Based on this result, the user can chose a grid cell discretization which provides an adequate level of connectivity.

The FracLips analysis of grid cell size addresses only one of a number of issues which need to be considered in selecting a grid cell size. Future research is needed to address these issues in grid specification:

• definition of permeability ellipse at grid block scale
• existence of EPM at grid block scale
• stochastic field parameters
• scale dependent permeability
• scale dependent dual porosity behaviors
• compartmentalization

2.2.6.1.2 Fracture System Porosity

The fracture system porosity f_F depends solely on the fracture network geometry and can be directly calculated as the product of the fracture intensity expressed as fracture area per unit volume (P₃₂) and the storage aperture of the fractures (e):

(Equation 2-23)

where:

f_F = fracture system porosi ty
V_F = fracture volume, [L³]
V_Cell = grid cell volume, [L³]
A_F = fracture area, [L²]
e = fracture storage aperture, [L]
P₃₂ = fracture intensity as area per unit volume, [L²/L³]

No information about the continuum model is necessary to calculate f_F, since it is independent of grid cell size as well as of main grid directions.

Because the fracture system porosity depends on the number of fractures per unit volume, the fracture size and the fracture aperture distribution, a different porosity has to be calculated for every area of the continuum model where these parameters differ. In FracLips, the fracture system porosity can be calculated for each grid cell, based on the fracturing in that cell. The primary issue in definition of fracture porosity from fracture intensity P₃₂ is the selection of an appropriate measure for storage aperture e. Possible measures include:

• aperture derived from transient hydraulic response
• mechanical aperture
• aperture derived from fracture permeability or transmissivity ("cubic law")
• aperture derived from geophysical measurements (gamma density, matrix porosity)
• correlations to fracture size and orientation.

The algorithm being implemented in FracLips assumes that fracture storage aperture, e, is assigned to fractures primarily based on transient hydraulic response, such that it reflect the storage capacity of the fracture system rather than the mechanical aperture or the ability of the fracture to carry flow.

2.2.6.1.3 Directional Fracture System Permeability

The permeability of the fracture system depends on the fracture intensity, the connectivity of the fracture network and the distribution of fracture transmissivities. An approach for approximate solution of the fracture system permeability has been developed by Oda (1984). This approach is illustrated in Figure 2-21. The approach of Oda is as follows (following Doolin and Mauldon, 1995):

Oda (1984) starts with the orientation of each fracture in a grid cell expressed as a unit normal vector n. Integrating the fractures over all of the unit normals N, Oda obtained a tensor N_ij describing the mass moment of inertia of fracture normals distributed over a unit sphere:

(Equation 2-24)

where n_i are the components (direction cosines) of a unit normal to the fracture n with respect to orthogonal reference axes x_i.

For a specific grid cell with know fracture areas A_k and transmissivities T_k, an empirical fracture tensor can be obtained by adding the individual fractures weighted by their area and transmissivity:

(Equation 2-25)

where V is the grid cell volume and k represents the k^th fracture of the N fractures in the grid cell.

Oda's permeability tensor is derived from F_ij by assuming that F_ij expresses fracture flow as a vector along the fracture's unit normal. Assuming that fractures are impermeable in a direction parallel to their unit normal, F_ij must be rotated into the planes of permeability

(Equation 2-26)

where F_kk is trF_ij.

The Oda (1984) solution has the advantage that it can be calculated without requiring flow simulations. However, it does not take fracture size and connectivity into account, and is therefore limited to well connected fracture networks. Two alternative approaches could be incorporated to FracLips are illustrated in Figure 2-22:

• pathways analysis: an approximate permeability tensor for a fracture network at a specified grid-cell scale can be obtained by carrying out a graph theory search through the network to identify the most permeable paths and to calculate their effective properties (Dershowitz et al., 1997)
• block scale flow simulation: an approximate permeability tensor for a fracture network at a specified grid-cell scale can be obtained by carrying out a series of flow simulations with boundary conditions which define unit gradients in each specified direction (see, e.g., Clemo and Smith, 1997; LaPointe et al., 1996). The key to this approach is the specification of appropriate boundary conditions. The flow solution can be obtained from DFN flow solvers such as MAFIC (Miller et al., 1997).

These approaches could be incorporated into future versions of FracLips.

2.2.6.1.4 Matrix-Fracture Interaction Parameters (Spacing and _-factor)

The typical fracture spacing describes the average distance between two fractures or expressed in a different way, the average thickness of a matrix block not disturbed by fractures¹. This parameter gives information about the average distance a fluid has to move before reaching the high permeability fracture system in which it is traveling to (or from) a well. Hence, it is a measure of the accessibility of the matrix system through the fracture system. Kazemi (1970) introduced the _-factor, which combines average distances in three perpendicular directions to describe this characteristic of fractured reservoirs.

The typical fracture spacing and the resulting _-factor are calculated in the exact directions of the continuum model grid to be used later. Hence, the main directions of the continuum model grid must be known at the time of this exercise. The fracture spacing is independent of the grid block size.

The average fracture spacing is related to various aspects of the fracture system (e.g., orientation, intensity). Therefore it is important to calculate these parameters for each grid cell. The calculation of fracture spacing and _-factor is illustrated in Figure 2-23. For more information, see Dershowitz et al. (1997). These calculations will be incorporated into FracLips based on algorithms developed for FraCluster.

2.2.6.1.5 Permeability Barriers

Some DP codes, including ECLipse, allow for the effects of heterogeneous connectivity through the definition of permeability barriers. Those barriers exist between cells which are not hydraulically connected to their neighboring cell (Figure 2-20). The use of permeability barriers facilitates the analysis of compartmentalized reservoirs, and provides some of the functionality of DFN model heterogeneous connectivity.

In FracLips, permeability barriers are identified by carrying out graph theory analysis of each grid cell to determine which neighboring cells are not hydraulically connected to this cell. In the ECLipse reservoir model the transmissibility between not connected cells is then set to zero (e.g., using the MULTX keyword).

The location of permeability barriers is dependent on the realization of the discrete fracture network, i.e., in every realization the location of those barriers will be different.

2.2.6.1.6 Active/Inactive Cell

Some DP codes, including ECLipse allow for the effects of heterogeneous connectivity through the definition of "inactive" cells. Inactive cells are defined as cells which are not hydraulically connected to any of the six neighboring cells (Figure 2-24). These cells play no role in the reservoir model except to serve as flow barriers. The use of inactive cells facilitates the analysis of compartmentalized reservoirs, and provides some of the functionality of DFN model heterogeneous connectivity.

In FracLips, inactive cells are generated by identified by carrying out graph theory analysis of each grid cell to determine which calls are not hydraulically connected to any neighboring cells. These cells can be marked as inactive for the ECLipse reservoir model.

2.2.6.1.7 Well Permeability Thickness Kh

The permeability thickness Kh (Figure 2-25) is one of the key parameter for the determination of well inflow, i.e., the `quality' of a well connection to the reservoir. For existing wells this parameter is derived from well test analysis.

In a porous media it describes the product of the completion interval thickness and the layer permeability. In fractured reservoirs it describes the sum of fracture transmissivities, i.e., the number of fractures intersecting the well and their transmissivity.

For predictions (i.e., new well locations) this parameter is not known and therefore needs to be calculated from the simulation grid data (i.e., completion data and layer permeability). The use of DFN models allows predicting this parameter also for fractured reservoirs. Furthermore this can be done for different well locations and well directions, thus allowing the modeler to optimize well placement.

For each new well connection the number of fracture intersections is computed by FraClips. The transmissivities of the intersecting fractures is then used to calculate the well Kh. Note that this parameter is derived strictly from geometry, i.e., no flow modeling is carried out. Moreover this parameter does not take into account the connectivity of the fracture network around the well (flow dimension). This can be done by calculating the wells productivity or injectivity index.

2.2.6.1.8 Well Productivity/Injectivity Index

The well productivity/injectivity index (Figure 2-26) is used to describe the connection between a well and the reservoir. The steady state index is defined as:

J = Q / Dp (Equation 2-27)

where:

Q =production/injection rate
Dp = pressure difference between well and reservoir.

For existing wells this parameter can be measured and it is often used directly in ECLipse to model well performance. For predictions (i.e., new well locations) this parameter is not known and therefore calculated from the existing grid data (e.g., permeability). The use of DFN models allows predicting this parameter for every location and well direction in the model. This is done by simulating steady state production or injection tests with MAFIC.

This approach involves the use of flow simulations and therefore is much more time consuming than calculating the well Kh. It does, however, take the fracture network connected to the well (flow dimension) into account.

2.2.7 Task 4.1.2: Well Testing Data Acquisition

During the quarter, Marathon collected and processed well test and hydraulic response data, and provided the data for posting on the WWW server. Well test and hydraulic responses provided by Marathon during the quarter are illustrated in Figures 2-27 and 2-28.

2.2.8 Task 5.1.2: WWW Site Updates

During the quarter, Golder Associates updated the World Wide Web (WWW) Server for the project. The update included the following postings:

• Fifth Quarterly Progress Report (Dershowitz et al., 1996)
• Additional data provided by Marathon (Figures 2-27 and 2-28)

2.2.9 Task 5.2.1 Progress Reports

The quarterly progress report for the period March 1, 1997 to June 30, 1997 was prepared during the quarter, and was posted to the WWW.

2.2.10 Task 5.2.2: Research Reports

No new research reports were completed during the quarter.

2.2.11 Task 5.2.3: Presentations

During the quarter, Bill Dershowitz attended the DOE Reservoir Research Conference in Houston, Texas, to present the paper, "Fractured Reservoir Discrete Feature Network Technologies." In addition, Bill Dershowitz, Tom Doe, and Paul LaPointe attended the US Rock Mechanics Symposium in New York, New York, to present the papers, "Flow Compartmentalization and Effective Permeability in 3D Fractal Fracture Networks", and "Analysis of Heterogeneously Connected Rock Masses by Forward Modeling of Fractional Dimension Flow Behavior."

2.2.12 Task 6: Management

On July 11, Bill Dershowitz met with Violetta Ivanova of MIT at MIT in Cambridge, Mass to provide technical review and supervision of ongoing MIT research. No other significant project management activities were carried out during the quarter.

1 Fracture in this context means fracture within the considered size range. The matrix itself may well include micro-fractures which often are responsible for reasonable matrix permeability.