- Copyright ©2011. The American Association of Petroleum Geologists. All rights reserved.
Structural elements of deformation-band fault zones are implemented as volumetrically expressed building blocks, that is, fault facies, in a series of synthetic reservoir geomodels and simulation models. The models are designed and built to reproduce a predefined range of fault system configuration, sedimentary facies configuration, and fault zone architecture. Using petrophysical properties derived from published field studies, the geomodel realizations are run in a reservoir simulator to monitor reservoir responses to variations in modeling factors. The modeled fault zones act as dual barrier-conduit systems, resulting in simulation models that can capture contrasting waterfront velocities, changes in waterfront geometries, and flow channelizing and bifurcation in the fault envelopes. The simulation models also show the development and sweep efficiency of bypassed oil and poorly swept regions because of the presence of the fault zones. Statistical analysis reveals that the fault facies modeling factors can be ranked according to impact on reservoir responses in the following descending order: fault core thickness, the type of displacement function, sedimentary facies configuration, the fraction of total fault throw accommodated by fault core and damage zones, fault system configuration, and maximum damage zone width. Fault core thickness is the most important factor because it governs the space available for fluid flow in the fault-dip direction. Other modeling factors affect the reservoir responses by controlling the geometry and continuity of fluid flow paths in the modeled fault zones.
2nd revised manuscript received April 25, 2010
Muhammad Fachri earned a B.Sc. degree in geophysics and an M.Sc. degree in geology from the Institute Technology Bandung (Indonesia). In 2006, he entered the University of Bergen, where he is currently a Ph.D. candidate. His scientific interests include faults and fluid flow, reservoir modeling, upgridding, and upscaling. He is currently employed as staff geologist at Weatherford Petroleum Consultants AS.
Jan Tveranger was awarded his D.Sc. degree in geology from the University of Bergen, Norway in 1995. He is currently employed as research manager of geosciences at the Centre for Integrated Petroleum Research (CIPR) in Bergen. His scientific interests include sedimentology, structural geology, and reservoir modeling.
Nestor Cardozo received his bachelor's degree in geology from the Universidad Nacional de Colombia and his Ph.D. in geology from Cornell University. He is an associate professor at the University of Stavanger, Department of Petroleum Engineering. His scientific interests cover faults, their related deformations, and their implementation in reservoir models.
Øystein Pettersen is a principal researcher and professor in computational mathematics (reservoir engineering and rock mechanics) at the Centre for Integrated Petroleum Research (CIPR) in Bergen. His main research areas are within rock mechanics modeling and simulation, and reservoir simulation. His former career includes 18 yr with Statoil in different positions, such as advising and specialist reservoir engineer, and leader of Statoil simulation network.
In the last few years, the need for a tool supplementing the conventional use of fault transmissibility multipliers (Manzocchi et al., 1999) has been argued by several workers (Flodin et al., 2001; Rivenæs and Dart, 2002; Tveranger et al., 2005; Berg and Øian, 2007; Fredman et al., 2007; Manzocchi et al., 2008). To fill this need, the fault facies modeling technique was conceived as a way to allow explicit three-dimensional (3-D) fault-envelope structures to be included in industrial reservoir models (Tveranger et al., 2004, 2005; Braathen et al., 2009). The procedure generates explicit fault envelope corner-point grids from conventional reservoir models (Syversveen et al., 2006). Subsequent steps populate this grid with realistic volumetrically expressed fault facies (i.e., geologic elements originating from tectonic deformation of the host rock). Classification schemes and measurements of fault facies derived from field observations (Braathen et al., 2009) are used during this process, either for defining large-scale fault facies (Syversveen et al., 2006; Skorstad et al., 2007; Soleng et al., 2007) or for defining outcrop-scale fault facies (Fredman et al., 2008). After petrophysical modeling and optional upscaling, the fault envelope grid is merged with the conventional model before reservoir simulation.
The adaptation of existing methods for stochastic modeling of sedimentary facies to fault facies modeling was demonstrated by Syversveen et al. (2006), Skorstad et al. (2007), Soleng et al. (2007), and Fredman et al. (2008). Fredman et al. (2008) document the method's capability of qualitatively reproducing outcrop-scale fault observations. Recent developments of the fault facies modeling technique are directed toward improvements in (1) generating structure-imitating models that reproduce outcrop-scale statistical measures of fault zones (partly shown in this article); (2) incorporating smaller structural elements, such as slip zones and slip surfaces, into fault facies models; and (3) upscaling of fault facies geomodels.
An important step when developing a new 3-D stochastic modeling technique, such as fault facies, is to conduct sensitivity and screening studies to investigate the influence of modeling factors on model outcomes and responses (cf. Jones et al., 1995; Seifert and Jensen, 1999; Fredman et al., 2007; Tveranger et al., 2008). By using an experimental design, this article investigates the impact of fault facies modeling constraints and parameters on geomodel realizations, modeled fluid flow behaviors, and reservoir responses. The modeling constraints include fault system and sedimentary facies configurations, which are implemented using synthetic models. To simplify the experimental design, we focus on fault zone geometric properties (e.g., fault core thickness, damage zone width, etc.) and use simplified bulk porosity and permeability values. The results of published field and subsurface studies are used to constrain the tested modeling parameters, porosity, and permeability values.
DEFORMATION-BAND FAULTS AND FAULT FACIES DESCRIPTION
Deformation-band faults (Aydin and Johnson, 1978) typically consist of a fault core accommodating most of the displacement and surrounding damage zones. In this article, we define the core of deformation-band faults as the region characterized by a high density of slip zones, shale smears (if shale beds are present), host rock lenses, and fault rock pockets (Foxford et al., 1998; Braathen et al., 2009). The hanging-wall and footwall damage zones are defined as surrounding the fault core and are characterized by relatively minor and more distributed deformation. Deformation bands are the major structural components of the damage zones. This subdivision follows, for example, Clausen et al. (2003), Berg and Skar (2005), and Flodin et al. (2005).
In this study, the structural elements of deformation-band faults are not represented as individual discrete fault facies. Instead, we represent the basic structural elements in the model as fault facies associations (Braathen et al., 2009). For example, in the fault core, a deformed sandstone lens can be viewed as a fault facies association consisting of deformation-band facies and sandstone facies bounded by slip zones. Lenses exhibiting different parent rocks and internal deformation intensity can be classified as separate fault facies associations (see figure 7 of Braathen et al., 2009).
Our modeling study uses a faulted sandstone-mudstone succession. We define four fault facies associations in the fault core, all originating from deformation of the host rock: (1) low-strained sandstone lenses (Sst-FC-L), (2) low-strained mudstone lenses (Sh-FC-L), (3) high-strained sandstone lenses (Sst-FC-H), and (4) high-strained mudstone lenses (Sh-FC-H). Detailed description of each facies association in the fault core is given in Table 1.
A similar approach is used for defining fault facies associations in the damage zone. Sandstone containing deformation-band clusters is viewed as a fault facies association consisting of deformation-band facies and sandstone facies and can be characterized by combining properties such as deformation-band density and the appearance and characteristics of slip surfaces. This classification has previously been used for characterizing the Big Hole fault damage zone in Utah that cuts Jurassic Navajo Sandstone (figures 3–5 of Shipton and Cowie, 2001). In this study, fault facies associations originating from sandstone in the damage zone are defined based on deformation-band frequency. Deformation-band frequency data from the Bartlett fault in Utah (Berg and Skar, 2005) are used as an analog. The fault cuts Jurassic Slick Rock Member of the Entrada Sandstone and the Jurassic Moab Member of the Curtis Formation. Deformation-band frequency in Berg and Skar (2005) is given as bands per meter, which could be easily adapted to fault-perpendicular grid cell dimensions in our synthetic model. We chose the frequency of 20 deformation bands per meter as the value for subdividing low-strained sandstone and high-strained sandstone (Table 1) because in the Moab Member, this value corresponds to the boundary between outer and inner damage zones. For subdividing unstrained sandstone and low-strained sandstone (Table 1), we chose the frequency of 1 deformation band per meter. Furthermore, we assume that the mudstone beds in the damage zone experienced ductile deformation and that the intensity of this ductile deformation is related to the deformation intensity in the adjacent sandstone beds. Thus, the following six fault facies associations can be defined in the damage zone and surrounding host rocks (Table 1): (1) unstrained sandstone (Sst-U), (2) unstrained mudstone (Mst-U), (3) low-strained sandstone (Sst-DZ-L), (4) low-strained mudstone (Mst-DZ-L), (5) high-strained sandstone (Sst-DZ-H), and (6) high-strained mudstone (Mst-DZ-H).
A reservoir zone is typically underlain and overlain by thick shale-dominated zones that form cap rock and floor rock seals (Sorkhabi and Tsuji, 2005), and thus the top and bottom of the reservoir zone are assumed to be no-flow boundaries during reservoir simulation. In the fault zones, this assumption is considered valid for production time scale if the cap and floor rocks deformation had occurred completely in a ductile manner. However, cap and floor rocks in the fault zones could deform in several different mechanisms, for example, ductile folding and local drag (Hesthammer and Fossen, 1998), shale smearing (Færseth, 2006); and brittle fracturing, faulting, and shale gouge formation (Braathen et al., 2009). These have significant impact on hydrocarbon migration and entrapment (Karlsen et al., 2004). Therefore, in the fault zone, the validity of no-flow boundary conditions in the cap and floor rocks depends on the prevailing faulting mechanism.
In fault facies modeling, gridding of the fault envelope captures the rock volume occupied partly by these deformed cap and floor rocks, and thus additional fault facies, that is, lower and upper, should be defined (Table 1). In this study, we assume lower and upper fault facies deformed by ductile folding and shale smearing so that they form continuous deformation domain boundaries with the reservoir zone. A special field analog of this situation can be found in the Moab fault (in Utah) and its segments, which are documented in detail by Davatzes and Aydin (2005).
This section presents the workflow used for generating the heuristic geomodels (fault facies and surrounding sedimentary facies models) used in screening studies. The models were created by an iterative process designed to capture a predefined range of sedimentary and structural heterogeneities, forming a matrix explained in the last subsection.
Basic Conventional Model Configurations and Fault Envelope Grids
The basic model measures 1000 × 1000 m (3280 × 3280 ft) with a thickness of 30 m (98 ft). A structured corner-point grid with a uniform cell resolution of 25 × 25 × 2 m (82 × 82 × 7 ft) (x, y, z) was used. Two fault system configurations were introduced; all faults are aligned to the grid (parallel to the y axis) and dip 60°:
One normal fault with displacement of 50 m (164 ft) (throw = 43 m [141 ft]). This configuration is labeled “I” in geomodel code (see the last subsection).
Two parallel normal faults, each with displacement of 25 m (82 ft) (throw = 21.5 m [71 ft]). This configuration is labeled “II.”
Two sedimentary facies configurations were used, both with a volumetric 60/40 split between sandstone and mudstone:
Five horizontal layers consisting of three sandstone beds and two mudstone beds. This configuration is labeled “5.”
Seven horizontal layers consisting of four sandstone beds and three mudstone beds. This configuration is labeled “7.”
For investigating the influence of fault facies modeling constraints on model outcomes and responses, four conventional reservoir models were created by combining the two sedimentary facies configurations with the two fault system configurations (Figure 1).
A I-5 model is used to illustrate the basic workflow (Figure 2). A fault envelope grid (Figure 2B) was created for each fault in the four conventional reservoir models (Figure 1) using the procedure described in Syversveen et al. (2006). The fault envelope width was chosen based on the expected maximum damage zone width. Fault throw (FT) was used to define the maximum damage zone width (hence, the width of the fault envelope) because both parameters are positively correlated (Beach et al., 1999; Shipton and Cowie, 2001). In the single-fault models (I-5 and I-7), the fault envelope was set to extend two cells into the footwall and two cells into the hanging wall of the fault, yielding a fault envelope grid 100 m (328 ft) wide. The resolution inside this envelope was refined by a factor of 25 × 3 × 3 (82 × 82 × 10 ft) (x, y, z). In the double-fault models (II-5 and II-7), the fault envelope width of each fault was set to 50 m (164 ft) (one cell width into both footwall and hanging wall of the fault), and the refinement factor is 25 × 3 × 2 (82 × 10 × 7 ft) (x, y, z). With these refinements, fault envelope grid cells in all geomodels have a dimension of approximately 1 × 8.3 × 2 m (3 × 27 × 7 ft) (x, y, z).
The size and aspect ratio of the cells was chosen based on the following considerations. First, many outcrop observations of deformation-band density in fault damage zones (Davatzes and Aydin, 2003; Berg and Skar, 2005) are recorded along scan lines using the number of deformation bands per meter as a unit of measure, hence, the fault facies models with a perpendicular-to-fault resolution of 1 m (3 ft) can easily be compared with the available field data. Second, the cell thickness in the conventional reservoir models (Figure 2) is 2 m (7 ft). To maintain the connectivity of the flow unit in the regions of small displacement fields (i.e., damage zones), the vertical grid refinement of the fault envelope was chosen to yield a vertical resolution of approximately 2 m (7 ft) or less. Third, in the deformation-band faults, the shortest correlation length can be expected to be found in the direction perpendicular to fault strike (here represented by the x direction). Thus, the cell length in the x direction was set smaller than the cell length in the y and z directions, that is, less than 2 m (7 ft) (1 m [3 ft] in our case). The longest correlation length can be expected to be found in the direction parallel to the fault strike (here represented by the y direction). Hence, the cell length in the y direction was set larger than the cell length in the z direction, that is, more than 2 m (7 ft) (8.3 m [27 ft] in our case).
The chosen cell-aspect ratio also conforms to the observed bulk permeability structure in typical deformation-band faults. In the direction perpendicular to fault strike, where grid cells have the shortest length, bulk permeability in a deformation-band fault is typically lower than the bulk permeability in the directions parallel to fault (Antonellini and Aydin, 1994; Caine et al., 1996, Sternlof et al., 2004; Sorkhabi and Tsuji, 2005). In the directions parallel to fault, the effects of structural connectivity (e.g., because of deformation-band interaction, crosscutting of synthetic and antithetic structures) in decreasing bulk permeability are assumed to be similar, and the bulk permeability is typically much higher than that in the perpendicular direction. However, the need to capture sedimentary flow units as previously discussed hampers the use of cells with similar lengths in fault-strike and fault-dip directions. Hence, only in the fault-strike direction can the length of cells be set much longer than the cell length in the fault-perpendicular direction. With the cell-aspect ratio previously discussed, the grid is well suited for reservoir simulation using bulk permeability values derived from field studies.
Setting Up the Displacement Function
The displacement function in the fault envelope was applied after (1) resampling sedimentary facies in the conventional grid into the fault envelope grid (Figure 2C) and (2) restoring the sedimentary facies distribution inside the fault envelope grid to its prefault configuration (Figure 2D). Displacement in deformation-band faults can be accommodated by several mechanisms. In this study, we assumed that the faults were developed entirely by shear displacement of deformation bands and slip surfaces, and thus the displacement function was parameterized based on this assumption. Different fault envelope structures in the geomodels were generated by varying displacement function variables. These variables include (Figure 3): (1) the maximum damage zone width as a function of FT; (2) the fault core thickness as a function of FT; (3) the fraction of total FT accommodated by the fault core (i.e., fault core throw percentage [FCTP]) and damage zones as functions of FT; and (4) the type of displacement function in damage zones and fault core (Figure 3B).
To investigate the influence of these parameters on model outcomes and responses, we chose two deterministic values or functions for each displacement function variable.
Variable 1, maximum damage zone width as a function of FT:
Maximum damage zone width proportional to 2× FT. This corresponds to a fault envelope width of 86 m (282 ft) in the single-fault models and 43 m (141 ft) for each fault in the double-fault models.
Maximum damage zone width proportional to 1.77× FT. This corresponds to a fault envelope width of 76 m (249 ft) in the single-fault models and 38 m (125 ft) for each fault in the double-fault models.
The chosen values for variable 1 are in the middle range of the damage zone width data set reported by Beach et al. (1999).
Variable 2, fault core thickness as a function of FT:
Fault core thickness proportional to 0.24× FT. This corresponds to a fault core thickness of 10.4 m (34 ft) in the single-fault models and 5.2 m (17 ft) for each fault in the double-fault models.
Fault core thickness proportional to 0.16× FT. This corresponds to a fault core thickness of 6.9 m (23 ft) in the single-fault models and 3.45 m (11 ft) for each fault in the double-fault models.
The chosen values for variable 2 are in the middle to upper range of the fault core thickness data (i.e., relatively thick fault core) published by Walsh et al. (1998) and Sperrevik et al. (2002). Outcrop analogs of isolated fault segments with very thick fault cores are shown in Shipton and Cowie (2001) (see their figure 9, where the fault core thickness is proportional to 0.25× FT) and Davatzes and Aydin (2003) (see their figure 16, where the fault core thickness is proportional to 0.14× FT).
Variable 3, the fraction of total FT accommodated by fault core and damage zones:
An FCTP of 90% (damage zones accommodate 10% of total FT).
An FCTP of 95% (damage zones accommodate 5% of total FT).
The chosen values for variable 3 were constrained by field and subsurface data from Davatzes and Aydin (2003) (for 2.6 and 110 m [9 and 360 ft] FT, see their figures 15, 16), Hesthammer and Fossen (2001) (for 6 and 69 m [20 and 226 ft] FT, the data in well 34/10-A-15 and well 34/10-B-12), and in-house Center for Integrated Petroleum Research (CIPR) data set (for 100-m [328 ft] FT). These data show that FCTP increases with increasing FT. We used a logarithmic function for determining FCTP as a function of FT: FCTP = 5.6046 ln(FT) + 71.633. This function conforms to the model of deformation-band fault evolution (Aydin and Johnson, 1978) because it imposes positive FT values when FCTP is zero, that is, a situation corresponding to the early stage of deformation-band formation and critical deformation-band density before slip surface formation. With this function, a fault with a 43 m (141 ft) throw corresponds to an FCTP of 93%, whereas a fault with a 21.5 m (70 ft) throw corresponds to an FCTP of 89%. Thus, the chosen values (90% and 95%) cover the range of FCTP values calculated based on selected field data.
Variable 4, the type of displacement function in damage zones and fault core:
Quadratic function in the damage zone (displacements are proportional to the fourth power of distance) and cubic function in the fault core (displacements are proportional to the third power of distance). This combined function is shown as a dotted black curve in Figure 3B and denoted as Q2-C.
Quartic function in the damage zone and quadratic function in the fault core (displacements are proportional to the second power of distance). This combined function is illustrated by a solid gray curve in Figure 3B and denoted as Q4-Q2.
The type of displacement function in the damage zone can be constrained using the cumulative frequency of deformation bands measured along scan lines perpendicular to the fault. For example, in several scan lines along the Bartlett fault (Berg and Skar, 2005), the gradients of curves representing the cumulative frequencies tend to increase toward the fault core. However, the increase in curve gradient varied spatially. In one scan line, the curve gradient increases more gradually, whereas in other scan lines, they increase more abruptly. If each deformation band is assumed to host a constant shear displacement, these curves can be used to derive a displacement function in the damage zone. Based on simple calculations, we find that the quadratic and quartic functions are good for capturing the range of displacement functions in the damage zone. The use of quartic function reflects a more localized deformation toward the fault core and consequently a thinner inner damage zone, whereas the use of quadratic function corresponds to a more distributed deformation and a wider inner damage zone. In the fault core, the quadratic function corresponds to a fault core with a single main slip surface at the hanging-wall side (Braathen et al., 2009) and represents a more localized deformation. The cubic function, however, corresponds to a fault core bounded by two main slip surfaces at both the footwall and hanging-wall sides (Hesthammer and Fossen, 2001; Davatzes and Aydin, 2003) and represents a more distributed deformation. In fault facies and strain modeling studies, linear and quadratic functions have been used by Soleng et al. (2007), Cardozo et al. (2008), and Fredman et al. (2008). However, the bases for determining the functions are different from those used in this article.
A combination of values from each of the four variables (1–4) forms one distinct displacement function. Using two alternative values for each variable yields a total of 16 possible displacement functions.
Applying the Displacement Functions, Probability Functions, and the Sequential Indicator Simulation of Fault Facies
The information of damage zone width and fault core thickness contained in the displacement function was used to create a fault-architectural-element parameter (Figure 2E). The displacement function was used to displace facies in the restored or resampled sedimentary facies parameter (Figure 2D), and the result of this procedure is a lithology distribution parameter (Figure 2F). Shear strain (Figure 2G) was computed by calculating the change in the angle of lines that were perpendicular before deformation (e.g., Ramsay and Huber, 1983). The shear strain calculation was conducted using the displacement gradient and local fault dip values provided by the displacement function and grid data. These parameters were subsequently used to create probability functions of fault facies (Figure 2H). During this step, fault-architectural element and lithology distribution parameters were used to constrain the location of fault facies, whereas the shear-strain parameter was used to constrain the spatial distribution of deformation intensity.
In the damage zones, the subdivision of the fault facies is based on the deformation-band frequency in Berg and Skar (2005) (see their figure 6). The shear strain parameter was used for capturing the spatial distribution of these fault facies. For example, both unstrained host rock and low-strained fault facies occur in the transitional zone, whereas unstrained host rock is absent in the outer zone. Because the minimum deformation-band frequency in the outer zone is one deformation band per meter, a shear strain value corresponding to this frequency was used to limit the distribution of unstrained host rock. We have used the simplest case for calculating this strain value: A vertical deformation band hosting 2.5 mm displacement along a 1-m (3.3-ft) distance corresponds to a shear strain of 0.0025. Thus, in the damage zones, we set the probability of unstrained fault facies equal to zero if the shear strain exceeds this value. A similar procedure was applied to the spatial distribution of low-strained fault facies and high-strained fault facies using a shear strain value corresponding to the frequency of 20 deformation bands per meter (Figure 2H, right).
In this study, the calculation of probability functions for the fault core facies follows the procedure used by Fredman et al. (2008) (see their figures 10, 11). They showed that because the displacement field in the fault core is more complex than that in the damage zone, the probability of fault core facies originating from a particular host rock such as sandstone was given a value larger than zero (10% in this case) in the mudstone part of the lithology distribution parameter (i.e., not only in the sandstone part such as for the damage zone). Figure 2H (left) displays the probability function of low-strained sandstone lenses in the fault core.
Fault facies probability functions (Figure 2H) were used as 3-D trends in the subsequent sequential indicator simulation (SIS) (Deutsch and Journel, 1992; Seifert and Jensen, 1999). Because the correlation length of fault facies has partly been represented by the chosen grid-cell-aspect ratio, the variogram range of fault facies was set equal to the dimensions of the grid cell, that is, 1 × 8.3 × 2 m (3.3 × 27.2 × 6.6 ft) (x, y, z). This is not a typical use of SIS because the only modeling parameters that were honored during the modeling are the fault facies probability functions. However, as shown in the next section, this choice ensures the reproduction of the desired fault zone structure in the resulting models: (1) highly alternating fault facies in the transitional damage zone and (2) the presence of lenses with varying sizes and internal deformations in the fault core. Figure 2I shows one of the resulting fault facies model after SIS. Fault facies models serve as input to the subsequent petrophysical modeling before merging of the fault envelope grid with the basic model (Figure 2J).
Fault facies modeling constraints and parameters were changed in a systematic manner in the basic workflow above to generate a series of geomodels. By applying 16 displacement functions in the 4 conventional reservoir models (Figure 1), 64 geomodels were generated. These geomodels are organized in a matrix form (Table 2). Each geomodel code consists of numbers indicating its original conventional reservoir model (e.g., I-5, II-7) and a number representing the applied displacement function. In the matrix (Table 2), the displacement function variables of each geomodel can be traced to the header of the matrix. Ten SIS realizations were executed from each model, yielding a total of 640 realizations. Subsequent reservoir simulations were conducted on all realizations.
FAULT FACIES MODEL ASSESSMENT AND CHARACTERISTICS
Models Versus Outcrop
Figure 4 shows a comparison between one generated fault facies model (Figure 4A) and a fault zone exposure. The damage zone model consists of unaffected protolith, transitional zone (consisting of unstrained host rock and low-strained fault facies), outer zone (consisting of low-strained fault facies), and inner zone (consisting of high-strained fault facies) (Figure 4B, right). The boundary between the outer and inner zones is sharp, and the width of the inner zone is constant, as suggested by field data from Utah (Berg and Skar, 2005) (Figure 4B, left) and western Sinai (Fredman, 2007). Variations in damage zone width also appear to be captured by the model, as indicated by two scan lines in Figure 4B (right). This variation is smaller than the typical variation reported in field studies that can be on the order of 25 m (82 ft) (Shipton and Cowie, 2001).
Fault facies in the fault core form clusters that mimic lenses entrained in the fault core (Figure 4A). Figure 4C (right) displays the distribution of sandstone and mudstone lenses in the model that resembles the distribution of lenses in natural deformation-band fault cores (Figure 4C, left). Lens frequency decreases with increasing lens size. The reproduction of the actual shapes of the lenses is limited by the grid resolution used in the model, but the overall lens geometry is represented, that is, lenses exhibit an elongated geometry both in fault-dip (Figure 4C, right) and fault-strike (Figure 4A) directions. Lense positions are linked closely to the applied displacement function. In Figure 4, because the cubic displacement function in the fault core was used for generating the model, high-strained facies are more distributed at the footwall and hanging-wall side of the fault core, whereas low-strained facies are more distributed in the middle part of the fault core. In general, the comparison of model outcomes and outcrop observations confirms that the model setup reproduces realistic results in terms of fault facies distributions at the chosen scale of modeling.
Comparing Fault Facies Model Outcomes with Different Parameters Settings
The organization of the geomodels in a matrix form (Table 2) allows us to investigate and compare the geomodels systematically. A straightforward way of doing this is to compare fault facies volumetric proportions in different parts of the grid (Figure 5). As shown in Figure 5A, when all other settings are kept equal, the volumetric proportions of fault facies in the fault envelope grids for single-fault configurations are virtually identical with those of double-fault configurations (Figure 5A). Thus, in the present geomodels, the effect of changing fault system configuration from single fault to double faults is basically to distribute the fault facies proportionally around the faults in the double-fault models. The zonal distribution patterns are similar regardless of the number of faults. The proportionality occurs because we applied linear relationships between FT and maximum damage zone width/fault core thickness. If nonlinear relationships are applied, it is expected that the resulting fault facies volumetric distributions will not be proportional.
Changing the configuration of sedimentary facies from five to seven layers has little impact on the distribution of fault facies volumetric proportions. As shown in Figure 5B, the patterns for both configurations are similar. Figure 5C shows that using a narrower damage zone width and a thinner fault core yields a higher volumetric proportion of high-strained fault facies. The effect of this modification is more pronounced in the fault core (Figure 5C, highlighted red bars). Increasing the displacement percentage in the fault core or damage zone produces a higher volumetric proportion of high-strained facies. The effect of this modification is more pronounced in the damage zones (Figure 5D, highlighted orange bars).
Using a quadratic displacement function in the damage zone yields a higher volumetric proportion of high-strained facies, whereas using a quartic function produces a lower volumetric proportion of high-strained facies (Figure 5E, highlighted orange bars). In the fault core, using a cubic function produces accumulation of high-strained facies in the footwall and hanging-wall side; using a quadratic function causes an accumulation of high-strained facies in the hanging-wall side (Figure 5F, highlighted red bars).
For a particular structural model, the spatial distribution of fault facies is mainly influenced by the displacement function parameters, the shear strain, and the fault facies probability functions. The present investigation (Figure 5) is a systematic and quantitative examination of these effects, which are required for relating geology to the reservoir simulation results of the geomodels discussed in the following sections.
Reservoir simulations are used to investigate the effect of fault facies modeling factors (i.e., fault system configuration, sedimentary facies configuration, and displacement function variables) on reservoir responses. For this purpose, we have assigned simple petrophysical properties to the sedimentary and fault facies and used a straightforward setup for production strategy. Fault zone flow behaviors are explained before the statistical analysis discussed in the next section.
Petrophysical Properties and Reservoir Simulation Setup
Constant bulk porosity and permeability values were assigned to each fault facies. For fault facies originating from sandstone, we postulated a systematic decrease in these values with increasing strain (Table 3). The difference in porosity values is not big (as compared with the porosity of unstrained sandstone) because we have chosen lowly deformed sandstone lenses as the most deformed fault facies (Table 1). With increasing strain, fault-normal bulk permeability (permX) decreases faster than fault-parallel bulk permeability (permY and permZ). As suggested by Antonellini and Aydin (1994) and Sternlof et al. (2004), the difference between fault-normal and fault-parallel bulk permeability can be more than two orders of magnitudes depending on volume fractions of deformation bands, the difference between host rock and deformation-band permeability, and the connectivity of deformation bands. In assigning permeability values, we assume that at the modeled reservoir, depth-confining pressure holds most of the slip surfaces shut; hence, their highly conduit effects on low-pressure fluid flow during production can be disregarded.
Mudstone is assumed to deform in a ductile manner, and thus the fault facies originating from mudstone were given identical (low) porosity and permeability values as compared with those of unstrained mudstone. Furthermore, because lower and upper facies are assumed to deform by ductile folding and shale smearing (see previous section on fault facies description, Table 1), they are still considered impermeable. The porosity and permeability values assigned to the fault facies are listed in Table 3.
All 640 geomodel realizations were transferred to a reservoir simulator. The dynamic properties used for the simulation are summarized in Table 4. One vertical water injector and one vertical producer were placed in the central parts of the hanging-wall and footwall blocks, respectively (Figure 6A). The connection between the wells and the reservoir was defined in cells containing sandstone facies. The liquid production was kept constant at 1,500 Sm3/day. The injection was controlled by 100% voidage replacement, that is, the reservoir volume water injection rate was controlled so that it equaled the production voidage rate. Each simulation was run until the water cut in the production well exceeded 90%. The simulations were run in parallel using 10 central processing units. The simulation times of all 640 runs ranged between 2 and 10 hr.
Fluid Flow Characteristics
Figures 6B to F show snapshots of the reservoir simulation in one of the model II-5-1 realizations (see Table 2). This model is characterized by two parallel normal faults (fault 1 in the east and fault 2 in the west) each with throw of 21.5 m that subdivide a zone of five sedimentary layers (three sandstone beds and two mudstone beds) into three blocks (Figure 6A). For each fault, the maximum damage zone width is 43 m (141 ft), and the fault core thickness is 5.2 m. A quartic displacement function is used in the damage zone and a quadratic function in the fault core. Figures 6B to F illustrate the flow characteristics inside the fault zone using oil saturation through elapsed time. From the beginning of water injection in the hanging-wall block, the fastest movement of the waterfront occurs in the upper sandstone, whereas the slowest movement occurs in the lower sandstone (see arrows in Figure 6B). The steepest waterfront occurs in the upper sandstone, whereas in the lower sandstone, the waterfront is gently dipping (see marked region in the lower sandstone in Figure 6B). The contrasting velocity and geometry of the waterfront in the different sandstone beds are caused by the configuration of the fault zone, especially the fault core configuration. The displacement function applied to this model causes a thickening of the highly permeable facies adjacent to the upper sandstone, whereas the highly permeable facies adjacent to the lower sandstone are thin. This configuration facilitates fluid displacement in the upper sandstone and causes substantial gravitational segregation of fluids in the lower sandstone in the areas between the injector well and fault 1.
The waterfront gets steeper as it moves into the fault 1 envelope (see marked regions in the middle and upper sandstones in Figure 6B). This is caused by the permeability anisotropy of the fault zone, where sandstone facies permeability is decreased in all directions as compared with the host rock, but the decrease is higher perpendicular to the fault compared with other directions (Table 3). The change from gentle to steep waterfront involves vertical movement of water. This may cause temporary flow channeling and bypassing of oil.
After passing through the core of fault 1, the waterfront in the upper sandstone in the hanging-wall block advances into the middle and upper sandstones but not into the lower sandstone in the middle block (Figure 6B, C). The flow into the middle sandstone is facilitated by the presence of sand in the fault facies suite derived from the upper mudstone bed, where we assigned a 10% probability for the presence of sandstone facies (see Figure 2H left for an analog). If we had not assigned a probability for sandstone facies to occur in the mainly mudstone part of the lithology distribution parameter, the waterfront would only have entered the upper sandstone in the middle block (flow to the same sandstone bed).
The waterfront in the upper sandstone in the hanging-wall block cannot flow to the lower sandstone in the middle block (Figure 6B, C) (although these two sandstone beds are juxtaposed) because of the presence of two shale smears that act as flow barriers. This flow characteristic shows that although we have distributed sandstone facies within the mudstone part of the lithology distribution parameter using a 10% probability, the number is not high enough to allow the waterfront to flow between the two sandstone beds crossing the two shale smears.
When crossing fault 1, the waterfront in the upper sandstone moves faster in slab 4 than in slab 3 (marked by ellipses in Figure 6C). This is unexpected because both the injection and production wells are located in slab 3. The cause for this flow behavior is the configuration of fault facies inside the fault core. In slab 4, high-permeability facies exhibit a high vertical continuity, whereas the low-permeability facies are laterally discontinuous (vertical arrow in slab 4 in Figure 6A). This allows the waterfront to move fast vertically in the fault core. In slab 3, mudstone and sandstone lenses (both low- and high-strain types) display a high lateral continuity (lateral arrow in slab 3 in Figure 6A), funneling the flow laterally inside the fault core and barring vertical movement (cf. long arrows in fault 1 in Figure 6C, D).
The flow in fault 2 is characterized by an oblique propagation of the waterfront (oblique arrows in Figure 6D). The flow behavior previously discussed causes oil to be bypassed in many places, especially in the damage zone of the upper sandstone (see black ellipses in Figure 6C–E). Snapshots taken later during the simulation show that this bypassing is temporary. This is especially evident in areas close to the injection well and in a direction perpendicular to the fault from the injection well (compare black ellipses in slab 2 in Figure 6D, E), where the water saturation is seen to increase considerably 1 month after the initial bypassing of oil. However, it can be expected that these flow processes, to some extent, may force the oil to flow away from the production well. Furthermore, in problematic areas that are far away or not in a direction perpendicular to the fault from the injection well, the sweep efficiency is very poor (compare black ellipses in slab 1 in Figure 6C–E).
Figure 6F is a simulation snapshot after 33 months. Figure 6F shows that the poorest oil sweep occurs in the lower sandstone in the hanging-wall block. Interestingly, the sweep efficiency in the lower sandstone in the middle block is much better, although both sandstone beds are connected to the faults (faults 1 and 2) with the same configuration (indicated by ellipses). The fluid front in the lower sandstone experiences gravitational segregation (and thus it becomes gentler) as it moves from the injector well toward fault 1, causing a poor sweep inside the envelope of this fault. Conversely, the waterfront in the lower sandstone is steep as it exits the fault 1 envelope, and because the two fault zones are close together, it is steep when it enters the fault 2 envelope, thus facilitating better sweep.
Reservoir simulations in other model realizations also showed similar flow characteristics as previously discussed. In summary, important flow characteristics are as follows: (1) contrasting waterfront velocities in different flow units; (2) changes in waterfront geometries; (3) flow channelizing and bifurcation (mainly in the fault core) in lateral, vertical, and oblique directions; (4) the development and sweep efficiency of bypassed oil regions; and (5) poorly swept regions caused by the presence of the fault zones.
For each reservoir simulation run, six reservoir responses were monitored: (1) time for injection well bottomhole pressure to reach 500 bars (T.500b), (2) time to water breakthrough (T.WBT), (3) time to abandonment (T.Abnd), (4) water cut after 790 days (WC.790d), (5) oil recovery after 790 days (OR.790d), (6) oil recovery at the end of simulation (OR.End).
The simulation results for all reliable runs are shown in Figure 7. We have excluded simulation results with nonlinear convergence problems; hence, the number or realizations in each model is not always 10. Out of 640, the total number of reliable simulation runs is 568. These problems mainly occurred in models with seven layers and thin fault core, where the probability to render realizations with very poor pressure communication across the fault is high. However, in models with a very good pressure communication such as models 3 and 11 of I-5 configuration, the injection well bottomhole pressure did not reach 500 bars, hence, the data are not available. The above contesting problems can be avoided using a more complex injection-production scheme, but this was not implemented because it would render a simulation data set that is difficult to analyze. Apart from the elimination of a minor number of simulation data, these problems provide early indications that the modified fault facies modeling factors have a huge impact on reservoir responses.
In most of the calculated responses, the models show a large variation in mean values and spread of outcomes (Figure 7). Only in OR.End is this variation small. From these results, we can quantitatively estimate the effect of fault facies modeling factors on the calculated reservoir responses (discussed in the next section).
Factorial Design at Two Levels
The geomodel matrix (Table 2) can be rearranged into the model setup shown in the bottom of Figure 7 (note that we have replaced the codes by colors). In experimental statistics, the models setup and associated reservoir simulation results are basically a 26 factorial design at two levels because we have selected two values or functions for each of the six modeling factors and run reservoir simulations for all possible combinations. Therefore, we use an established method for analyzing factorial design at two levels (Box et al., 2005, ch. 5) to perform screening studies on fault facies modeling factors. This method involves the calculation of the “main effect” (of each factor) and the “interaction effect” (between two or more factors) on a particular response. We explain the definition and calculation of these effects with examples as follows.
Figure 8A (left) illustrates the calculation of the main effect of displacement function type on T.500b. The difference in average T.500b caused by localizing deformation (changing displacement function type from Q2-C to Q4-Q2) when the remaining factors are held fixed supplies one measure of the effect of displacement function type on T.500b (one point in Figure 8A, left). With 26 (64) models, 32 such measures are present. The average of these 32 measures is called the “main effect” of displacement function type on T.500b (indicated by a red line). Similarly, Figure 8A (right) illustrates the calculation of the main effect of FCTP on T.500b. This calculation procedure was conducted for all factors on all responses. Comparison between the two diagrams in Figure 8A shows that the effect on T500b caused by variations in displacement function type is larger than that because of variations in fault core thickness percentage.
In most cases, a main effect is not a complete description of the influence of a factor on a particular response. A factor can interact with two or more factors and have a combined influence upon a response beyond their main effects. Figure 8B (left) illustrates the calculation of a two-factor interaction between fault core thickness and sedimentary facies configuration on T.WBT. As shown in this figure, the effect of decreasing fault core thickness on T.WBT is smaller in five-layer models than that in seven-layer models (indicated by dashed blue lines). Hence, the factors fault core thickness and sedimentary facies configuration are said to interact. Half of the difference between the two effects is called the interaction effect (indicated by a solid red line). Similarly, Figure 8B (right) illustrates the calculation of a two-factor interaction between sedimentary facies configuration and fault system configuration on T.WBT. A comparison between the two diagrams in Figure 8B shows that the interaction effect of fault core thickness and sedimentary facies configuration on T.WBT is larger than that of sedimentary facies configuration and fault system configuration. Higher order interactions are calculated similarly. With six factors, we can have up to six factor interactions. This calculation procedure was conducted for all factor interactions on all responses.
To determine which effects are almost certainly real and which might readily be explained by stochastic variations, the standard error (SE) of effects needs to be estimated. By assuming that the ratio between the effect and the SE of the effect is normally, identically and independently distributed, within a 95% confidence interval, effects greater than 1.997 (≈ 2) times their SE are considered to be almost certainly real (real effects), that is, they are almost certainly not caused by stochastic variation. Table 5 shows some of the estimated effects along with their SEs. In Table 5, real effects are shown in bold type, whereas effects greater than one-tenth of the highest effect in each reservoir response (significant effects) are shaded gray. Because T.500b data for models 3 and 11 of the I-5 configuration were not available, we only calculated the main effects for this reservoir response.
Table 5 shows that, except for OR.End, all the other reservoir responses are related, especially for significant effects (gray shaded). Based on their trend, the responses can be subdivided into two groups: (1) T.500b, WC.790d, and OR.790d and (2) T.WBT and T.Abnd.
This response subdivision is used for subsequent statistical inferences as follows. For the main effects, the effects of displacement function, FCTP, fault core thickness, and sedimentary facies configuration decrease group 1 responses and increase group 2 responses, whereas the effects of fault system configuration increase group 1 responses and decrease group 2 responses. Except for OR.End, the main effects of displacement function and fault core thickness are significantly higher than those of the other factors. In general, displacement function is more important in responses related to a previous simulation time (T.500b and T.WBT), whereas fault core thickness is more important in responses related to later simulation time (T.Abnd, OR.790d, and OR.End).
The main effects of modeling factors cannot be quantified separately because of the factor interactions. Except for OR.End, significant effects of interactions occur up to three factors, whereas for more than three-factor interactions, the effects are small. In two-factor interactions, the interaction of displacement function and fault core thickness yields the greatest effect, whereas for three-factor interactions, the largest effect is shown by the interactions of displacement function, fault core thickness, and sedimentary facies configuration. The modification of maximum damage zone width yields the smallest effect. Because neither the main effect of maximum damage zone width nor any of the interactions involving maximum damage zone width (e.g., 1 × 4, 1 × 3 × 4, etc.) are real, maximum damage zone width is inert (at least for the factors and levels tested).
Table 5 shows that OR.End is not related to the other responses. Compared with the other responses, OR.End shows a higher number of significant interactions (numbers in bold type) that occur up to five factors (i.e., most of the effects are greater than one-tenth of the highest effect, −2.0 ± 0.1). The main effects of this response are considerably smaller than those of similar responses (i.e., OR.790d). This suggests that no factors significantly control OR.End.
The ranks in each response in Table 5 were given based on the absolute value of effects. Using (1) the ranks, (2) the number and magnitude of significant interactions involving the factors, and (3) the overall importance of factors across all the responses, the tested fault facies modeling factors (at least for the levels used) can be graded in the following descending order:
Fault core thickness
The type of displacement function
Sedimentary facies configuration
Fault system configuration
Maximum damage zone width
Effects of Internal Fault Zone Architecture
Because faults are represented as volumetric entities, the reservoir responses mainly depend on how easy and effective the fluids flow along the fault-dip direction. Fault core thickness is the most influential modeling parameter governing the flow characteristics because it affects the space available for flow in the fault-dip direction. This influence is strengthened by the fact that thinning of the fault core results in an increase in the volumetric proportion of high-strained low-permeability fault facies in the fault core (Figure 5C, red bars). The same judgment can, in other cases, be given to the maximum damage zone width parameter. However, because the present modeled geology (ductile mudstone deformation) leads to the implementation of strict vertical segmentation of flow units in the damage zones, the occurrence of vertical fluid flow within the damage zones is minimized, and consequently, the maximum damage zone width turns out to be the least important parameter affecting the reservoir responses.
As illustrated in Figure 6 and the previous discussion, the modeled fault cores impede fluid flow across the faults and form preferential paths for flow along the faults (in lateral, vertical, and oblique directions). Faults with such behaviors are said to act as dual conduit-barrier systems to fluid flow (Bense and Person, 2006). Faults in the Bighorn Basin in Wyoming (Bredehoeft et al., 1992) and the Baton Rouge fault in southern Louisiana (Stoessell and Prochaska, 2005) are among some of the faults in siliciclastic sedimentary rocks that show evidences of a conduit-barrier system. Bredehoeft et al. (1992) used history-matching techniques (by adjusting model parameters) and concluded that the volumetric representation of faults as a conduit-barrier system was needed to match the simulated and interpreted potentiometric map of the Bighorn Basin. Furthermore, hydraulic head, geochemical, and isotope data in the area covering the Baton Rouge fault provide more direct evidences that might be used to classify the fault as a conduit-barrier system. This fault is known as a major barrier to groundwater flow, separating freshwater aquifers north of the fault and saltwater aquifers south of the fault (Miller, 2005). However, the fault has been interpreted as a conduit for brackish groundwater that originated from deeper formations to encroach fresh groundwater in shallow aquifers in the area just north of the fault (Stoessell and Prochaska, 2005). Similar to our study, Bense and Person (2006) demonstrated that the dual barrier-conduit system of the Baton Rouge fault can be captured in fluid flow models by representing fault zones with discrete grid cells and a strongly anisotropic hydraulic structure.
The statistical analysis (Table 5) shows that among the displacement function variables, the displacement function type and FCTP are the most important factors after fault core thickness. The application of these modeling parameters displaces the original host rock facies in the fault envelope (thus affecting local host rock facies connections) (Figure 2F) and distributes fault facies according to shear strain (Figure 5D–F). Therefore, these parameters influence the reservoir responses by governing the geometry and continuity of flow path in directions perpendicular and parallel to fault dip (e.g., high-permeability lenses in the fault core). In summary, compared with their counterpart models, models with thinner fault core, models with Q4-Q2 displacement function type, and models with higher FCTP correspond to fault zones with more localized deformation. During the modeling, a more localized deformation corresponds to a more limited space for vertical flow and a tendency for creating noncontinuous flow path in directions perpendicular and parallel to fault dip. As a consequence, these models yield faster (lower) group 1 responses (T.500b, WC.790d, and OR.790d) and slower group 2 responses (T.WBT and T.Abnd).
Along a natural deformation-band fault, regions with relatively more distributed and more localized deformations occur alternately (Foxford et al., 1998; Shipton and Cowie, 2003). For example, in a fault segment with a similar total throw, the fault core can be very thick and contains several stacking host rock lenses (indicating a more distributed deformation) that change gradually to a single slip surface with a very thin (localized) zone of deformation (e.g., Shipton and Cowie, 2001). Using the flow analysis previously discussed (Figure 6, Table 5) as an insight, we can expect the tendency of fluid front to flow to regions with more distributed deformation and, as a consequence, that the fluid front in fault zones is likely to be irregular. This hypothetical relationship between complex fault zone fluid flow and overall deformation intensity might be used to explain the different location and irregular spatial distribution of fluid front between brackish and fresh groundwater in shallow aquifers in the area just north of the Baton Rouge fault (Stoessell and Prochaska, 2005). The complexity of fluid flow in the fault-strike direction, in addition to the groundwater well location, was not addressed in the Bense and Person (2006) study. Based on the reservoir simulation results in our study, it is expected that a 3-D reservoir simulation models using volumetric representation of fault zone structure (e.g., fault facies) should be capable of capturing 3-D complexity of fluid flow in fault zones.
Effects of Fault System and Sedimentary Facies Configurations
The effect of changing the fault system configuration from a single fault to double faults is to distribute the displacement and fault facies around the two faults (Figure 5A). Thus, the spatial distribution of fault facies is changed significantly, and a large impact on reservoir responses is expected. However, because this modification distributes the displacement and fault facies around the two faults proportionally (caused by the use of linear relationships), it does not affect the total volume of space available for vertical fluid flow along the fault cores. Furthermore, the geometry and continuity of flow path in directions perpendicular and parallel to fault dip are kept similar (they are only proportionally allocated). This might explain why the modification of fault system from single fault to double faults, although it significantly changes the fault facies architecture, has little effect on reservoir responses (see Table 5 and its interpretation). Field data points relating fault core thickness and FT (figure 2 of Sperrevik et al., 2002) suggest that the relationship follows a power-law relationship. If, instead of using a linear relationship, we use a power function, then each fault in the double-fault models will have a thicker fault core. In this case, the partitioning of deformation of a single fault into deformation related to two faults is similar to making the fault zone deformation less localized. As discussed in a previous subsection, this modification will give models that yield slower (higher) group 1 responses and faster group 2 responses.
In relation with the fault system configuration, we extend and generalize the discussion to the flow characteristics (Figure 6F). The sweep efficiency in the lowermost sandstone in the middle block is much better than in the lowermost sandstone in the hanging-wall block, although both sandstone beds are connected to the faults (eastern and western faults) with the same architecture. The western damage zone of the eastern fault zone is responsible for the better sweep in the middle block by steeping the waterfront before it enters the western fault zone. As a consequence, double-fault models yield slower (higher) group 1 responses and faster group 2 responses than those of the corresponding single-fault models. In general, this flow analysis illustrates the dependence of fault block sweep efficiency on damage zone width. However, our statistical analysis of the effect of maximum damage zone width on reservoir responses (Table 5) shows that a 10-m (33 ft) change of damage zone width in the middle block gave a subtle effect. With extremely narrow damage zones, significant fluid gravitational segregation occurs in the fault block, and the fluid front is very gently dipping when it enters the fault block. Hence, contrary to our reservoir simulation results, it is expected that a fault system containing one or several fault blocks but with extremely narrow damage zones will consistently yield faster (lower) group 1 responses and slower group 2 responses than those of the corresponding fault system without fault blocks (i.e., single-fault models).
The effects of sedimentary facies configurations are discussed as follows. The patterns of the fault facies of five-layer and seven-layer models are similar (Figure 5B). However, the statistical analysis (Table 5 and its interpretation) reveals that sedimentary facies configuration is the third important parameter after fault core thickness and displacement function type. The modification of sedimentary facies configuration corresponds to thinning of the sandstone beds (but increasing their number). In the present fault facies modeling, two main effects of thinning of the sandstone beds occur: (1) the average size of lenses in the fault core becomes smaller, and consequently (2), the number of lenses in the fault core becomes higher. Fredman et al. (2007) conducted a sensitivity study on the effect of fault core architecture and permeability on simulated fluid flow. The modified parameters are fault rock permeability, sand lens permeability, sand lens fraction, and sand lens connectivity. This study concluded that the most important parameters influencing reservoir performance were fault rock permeability and whether the sand lenses were connected to the unstrained host rock. The size distribution of lens was not considered in this study.
As previously discussed, in this study, the architecture of fault core depends on the sedimentary facies configuration. The change of sedimentary facies configuration from five layers to seven layers does not affect the sand lens proportion and the total volume of space available for vertical fluid flow along the fault cores. However, because sandstone lenses in the fault core become smaller, the continuity of flow path in directions perpendicular and parallel to fault dip is more limited. Consequently, models with seven layers yield faster (lower) group 1 responses (T.500b, WC.790d, and OR.790d) and slower group 2 responses (T.WBT and T.Abnd). As a complement to the study by Fredman et al. (2007), our observations show that the size distribution of permeable host rock lenses, which in our study depends on sedimentary facies configuration, is also another important parameter controlling the reservoir responses beside fault rock permeability and lense-damage zone connection.
Expanding the Geomodel Matrix
One of the advantages of using factorial design at two levels is the simplicity of implementing a more detailed, but commonly limited, investigation (i.e., sequential assembly) (Box et al., 2005). Using sequential assembly, further tests can be conducted on maximum damage zone width. This factor showed no real effects, but because the range of the levels chosen is narrow, we cannot conclude whether this factor is essentially inert. Apart from the factor of maximum damage zone width, within the framework of this sensitivity study, the influence of other factors can be easily explained (see the previous discussions). Using the framework of the present fault facies modeling, a sequential assembly also can be conducted to test other factors by modifying the levels of previous factors and applying additional constraints. This is beneficial because adding more factors for another full factorial design leads to a much larger matrix. For example, the effect of fault-related folding can be tested using a much lower FCTP and by setting the resulting fault facies the same as those in the associated models with 90% FCTP. In addition, the effect of soft-linked fault interaction in two-fault models can be tested using a wider damage zone in the area between the faults.
Because introducing too many factors for a full factorial design leads to a much larger matrix and, as a consequence, cumbersome statistical analysis, in this study, we focus our investigation on geometric properties of fault zones. Another modeling factor that needs to be fully investigated is fault zone permeability. For this purpose, a full factorial design developed from the previous design (Table 2) needs to be implemented. After this study, at least two important tests occurred regarding the influence of fault zone permeability in fault facies modeling. First, another set of permeability values, as opposed to those shown in Table 3, should be tested. Using more isotropic and smaller bulk permeability values will likely yield very different reservoir simulation results than those of this study. Although several studies concerning bulk permeability in fault zones (Antonellini and Aydin, 1994; Sorkhabi and Tsuji, 2005) show highly anisotropic nature of fault zone bulk permeability, several faulting mechanisms can result in a more isotropic permeability structure. These mechanisms include intensive formation of high-angle antithetic structures and bed rotation of host rock lenses in fault core.
Second, in our study, the permeability of slip zones are represented by averaging them with the permeability of adjacent lenses. This procedure, which can be called permeability modification, would probably yield inaccurate simulation results. Therefore, a test on the slip zone representation should be conducted. Besides the permeability modification procedures used in our study, two different procedures that follow fault facies modeling can be used for capturing the effect of slip zones to fluid flow: (1) transmissibility modification between grid cells separated by slip zones and (2) explicit representation using high-resolution grid cells created using local grid refinement. The latter offers possibilities for detailed studies targeting the impact of highly permeable slip surfaces on fluid flow.
Upgridding and Upscaling of Fault Facies Geomodels
A major concern with respect to using the fault facies method for industrial purposes is the apparent need to add high-resolution fault envelope grids to reservoir models (Rivenæs and Dart, 2002; Manzocchi et al., 2008). Adding a substantial number of cells to a model comes at a high cost in terms of central processing unit time, which in extreme cases can render models unmanageable. Commonly, this issue is handled by upscaling high-resolution models into coarser grids that reduce the number of cells and replace detailed heterogeneities distributed over several cells in the original model with averaged values in larger cells occupying the same volume. Although being a well-established technique for handling high-resolution sedimentary facies models, the extent of fine-scale 3-D heterogeneity in fault zones offers additional challenges to upscaling. Initial tests using conceptual fault facies models (Soleng et al., 2007) suggest that fluid flow patterns observed in high-resolution models are reproduced during upscaling. However, this needs to be corroborated using realistic fault zone models and comparing the response of high-resolution models with upscaled versions to investigate the response of different simulation parameters to upscaling and establish practical guidelines.
This issue is addressed in an ongoing study within the framework of the Fault Facies Project, which proposes a combined upgridding-upscaling procedure for generating accurate coarse models derived from fine-scale fault facies models. The fine-scale models in this study are similar to the models presented in this study. Adapting the results of upgridding studies in the realm of sedimentary facies (Garcia et al., 1992; Durlofsky et al., 1997; Wen and Gomez-Hernandez, 1998; He and Durlofsky, 2006), the coarse grid is generated using some information from the fine-scale fault facies models. Using the grid, the assessment was conducted toward the effects of different coarsening factors and local boundary conditions on the accuracy of the generated coarse models. To evaluate their robustness, the upscaling procedures are tested against different fault facies modeling constraints: (1) the number of sedimentary layers and (2) FT. The experiment results show that coarse models using open boundary conditions yield satisfactory accuracy, although the coarsening factor is very high (2–3 orders of magnitude). In addition, the directions for reservoir simulation gridding are discussed. The results of this study provide a foundation for the generation of fieldscale reservoir simulation models from fault facies geomodels such as used in the present article.
As part of the development of a new volumetric fault modeling technique (fault facies), this article addressed the screening studies on the effect of modeling factors (fault system configuration, sedimentary facies configuration, and displacement function variables) to the resulting fault envelope structure and reservoir responses. Major conclusions of this study can be drawn as follows:
Modifying fault system configuration has little effect on the total volumetric proportion of fault facies if linear relationships between FT and fault core thickness or maximum damage zone width are used. Modifying sedimentary facies configuration also yields a small effect on the fault facies volumetric proportion. However, the modification of sedimentary facies configuration by thinning sandstone beds decreases the average size of lenses in the fault core and, consequently, the number of lenses becomes higher.
A narrower damage zone, thinner fault core, and higher displacement percentage correspond to a higher volumetric proportion of high-strained fault facies.
The modeled fault zones act as dual barrier-conduit system, and it is demonstrated that the simulation models can capture contrasting waterfront velocities, changes in waterfront geometries, and flow channelizing and bifurcation in the fault envelopes. The simulation models also show the development and sweep efficiency of bypassed oil regions, and poorly swept regions caused by the presence of the fault zones.
Fault facies modeling factors can be graded according to their impact on modeling outcomes and reservoir responses in the following descending order: fault core thickness, the type of displacement function, sedimentary facies configuration, FCTP, fault system configuration, and maximum damage zone width.
In the simulated reservoirs, the reservoir responses depend mainly on the space available for flow in the fault-dip direction and this is governed solely by the fault core thickness. The reservoir responses also depend on the geometry and continuity of flow paths in the modeled fault zone. These, to a great extent, are controlled by sedimentary facies configuration and, to a lesser extent, by the type of displacement function, FCTP, fault system configuration, and maximum damage zone width.
Within the framework of the present fault facies modeling, sequential assembly can be conducted to test other levels and factors such as fault-related folding and fault interaction. This procedure works by modifying the levels of existing factors and applying additional modeling constraints. Further tests on the effect of permeability values need to address two issues: petrophysical values and the representation of slip zones.
This research was conducted as part of the Fault Facies Project at the Centre for Integrated Petroleum Research (CIPR), University of Bergen, Norway. The Research Council of Norway, CIPR, Statoil, and ConocoPhillips provided financial support for this work. We thank David Moreno for discussion on the use of parallel computing in reservoir simulation. Silje S. Berg is thanked for permission to reproduce sketches in Figure 4C from her Ph.D. dissertation. The manuscript has benefited greatly from critical reviews by Eric A. Flodin, Rosalind A. Archer, and an anonymous referee. Roxar, the Norwegian Computing Center and Schlumberger generously provided the following software: Irap_RMS (for geomodeling platform), HAVANA (for fault envelope gridding) and Eclipse 100 (for reservoir simulation).
The AAPG Editor thanks the following reviewers for their work on this paper: Rosalind A. Archer, Eric A. Flodin, and an anonymous reviewer.
- Manuscript receivedAugust 3, 2009.
- Revised manuscript receivedOctober 9, 2009.
- Revised manuscript receivedJanuary 20, 2010.
- Final acceptanceSeptember 13, 2010.