Influence of deformation-band fault damage zone on reservoir performance

Access to 3D descriptions of fault zone architectures and recent development of modeling techniques allowing explicit rendering of these features in reservoir models, provide a new tool for detailed implementation of fault zone properties. Our aim is to assess how explicit rendering of fault zone architecture and properties affects performance of fluid flow simulation models. The test models use a fault with a maximum 100 m displacement and a fault damage zone with petrophysical heterogeneity caused by the presence of deformation bands. The distribution pattern of deformation bands in fault damage zones is well-documented, which allows generation of realistic models. A multiscale modeling workflow is applied to incorporate these features into reservoir models. Model input parameters were modulated to provide a range of property distributions, and the interplay between the modeling parameters and reservoir performance was analyzed. The influence of deformation-band damage zone on reservoir performance in the presence of different fault core transmissibility-multipliers was investigated. Two configurations are considered: one in which the fault terminates inside themodel domain, representing a case in which the fluid can flow around the fault, and one in which the fault dissects the entire model domain, representing a case in which the fluid is forced to cross the fault. We observed that the impact of deformation-band fault damage zone on reservoir performance changes when the fault core transmissibility multiplier is changed. Reservoir performance is insensitive to changing damage zone heterogeneity in a configuration in which flow can move around the fault. Where flow cannot bypass the fault, the influence of fault damage zone heterogeneity on reservoir performance is significant even when the fault core transmissibility multiplier is low. Introduction Since the 1970s, fault damage zones dominated by millimeterto centimeter-thick deformation bands have been identified in outcrops and core data worldwide (Aydin, 1978; Antonellini and Aydin, 1994, 1995; Fossen and Hesthammer, 1997; Shipton and Cowie, 2001, 2003; Fossen and Bale, 2007; Fossen et al., 2007; Cilona et al., 2012; Schueller et al., 2013). A significant number of studies have been dedicated to understanding thickness, length, displacement, orientation, and petrophysical properties of deformation bands (for review, see Fossen et al., 2007; Fossen and Bale, 2007). Other workers have addressed the scaling relationship between damage zone width and fault displacement (Knott et al., 1996; Beach et al., 1999; Fossen and Hesthammer, 2000; Shipton and Cowie, 2001; Schueller et al., 2013), and how deformation band properties are controlled by tectonic regime, host rock lithology, and burial depth (Ballas et al., 2014, 2015; Soliva et al., 2016). Deformation bands are characterized by features, such as pore-space collapse, grain rearrangement, and grain-crushing, which affect porosity and permeability, and can consequently, depending on the dominant mechanism of deformation, be expected to influence fluid flow in reservoirs by forming barriers, baffles, or conduits (Aydin et al., 2006; Fossen et al., 2007). The impact of individual deformation bands on fluid flow in sandstones has been extensively studied using physical experiments on core plugs and samples. A review of these studies by Fossen and Bale (2007) shows that, deformation bands, either individually or as clusters, may display a permeability reduction relative to their host rock of up to five to six orders of magnitude. However, the same study concludes that the impact on effective permeability is limited, and controlled by the cumulative thickness of deformation bands and their 3D connectivity. Understanding the influence of deformation-band fault zones on a reservoir scale requires proper implementaUni Research CIPR, Bergen, Norway and University of Bergen, Department of Earth Science, Bergen, Norway. E-mail: qudongfang2012@gmail .com; jan.tveranger@uni.no. Qarat Al Milh Petroleum, Muscat, Sultanate of Oman. E-mail: muhammad.fachri@qampetroleum.com. Manuscript received by the Editor 4 December 2016; revised manuscript received 31 March 2017; published ahead of production 18 May 2017; published online 14 July 2017. This paper appears in Interpretation, Vol. 5, No. 4 (November 2017); p. SP41–SP56, 13 FIGS., 4 TABLES. http://dx.doi.org/10.1190/INT-2016-0229.1. © The Authors.Published by the Society of Exploration Geophysicists and the American Association of Petroleum Geologists. All article content, except where otherwise noted (including republished material), is licensed under a Creative Commons Attribution 4.0 Unported License (CC BY). See http://creativecommons.org/licenses/by/4.0/. Distribution or reproduction of this work in whole or in part commercially or noncommercially requires full attribution of the original publication, including its digital object identifier (DOI). t Special section: Fault damage zones Interpretation / November 2017 SP41 Downloaded from https://pubs.geoscienceworld.org/interpretation/article-pdf/5/4/SP41/4059949/int-2016-0229.1.pdf by guest on 23 July 2019 tion of these features in reservoir models. In general, fault zones are conceived as consisting of a fault core accommodating the bulk displacement and a surrounding damage zone (Caine et al., 1996; Shipton et al., 2005). Conventional reservoir models use transmissibility multipliers to capture the effect of faults zones on cross-fault fluid flow. These are based on the calculated fault zone permeability, fault zone thickness, and permeability of the footwall and hanging-wall grid cells. Fault zone thickness is estimated from correlations against fault throw. The calculation of fault zone permeability is based on an empirical relationship between shale content and displacement (Manzocchi et al., 1999). The “fault zone thickness” in this context is defined as either the separation between the outermost slip surfaces minus the thickness of undeformed lenses, or the thickness of the slip-surface itself in a lacuna (Manzocchi et al., 1999). Identifying the “outmost slip surfaces” is an ambiguous exercise even in well-exposed outcrops. In practice, fault zone thickness is handled as an approximation of the fault core thickness and envisaged as a homogeneous volume of rock representing the deformed rocks present in the fault zone. It is difficult to make a case for also including the damage zone in fault zone thickness because damage zones display clear spatial trends in deformation intensity and commonly are anything but homogenized. Damage zones are up to several orders of magnitude wider than the fault core and extend tens of meters away from the core (Shipton et al., 2002, 2005; Torabi and Berg, 2011). For common reservoir models using 10–100 m XY resolution, damage zone features, if implemented, can be expected to influence fluid flow patterns. Some researchers captured the effect of damage zones by reducing the transmissibility multipliers of the nearby fault connections (Manzocchi et al., 2008). However, the resulting reservoir simulation models do not render explicit fault damage zone heterogeneity and accurate fluid flow calculations in the damage zones. Several studies have attempted to capture the effect of deformation-band damage zones on fluid flow using reservoir models with explicit representation of damage zones. Some of these studies focus on deriving bulk permeability values for implementation in low-resolution simulation models (Odling et al., 2004; Sternlof et al., 2004; Harris et al., 2007; Fachri et al., 2013b); hence, the spatial scale of their models is generally limited. Similar studies with larger model dimension (Rotevatn et al., 2009; Fachri et al., 2013a) are based on field observation in deriving deterministic upscaled damage zone permeability, leaving room for improvement by using stochastic method where input (e.g., deformation band orientation and volumetric proportion) can be modulated. It has not been clear whether, or under which circumstances, the spatial variation of fault zone properties influences regional flow patterns and reservoir behavior or not. The objective of this paper is to (1) demonstrate an improved workflow in damage zone modeling and (2) investigate the contribution of fault damage zone on reservoir performance and the sensitivity of reservoir performance to fault damage zone heterogeneity, thereby to infer implications regarding fault zone property modeling issues in real-field geologic models. For this study, we use a multiscale modeling workflow based around the use of fault facies (Tveranger et al., 2005; Braathen et al., 2009) to model a conceptual normal fault in porous sandstone with a damage zone consisting of deformation bands (for the previous applications of fault facies, see Fredman et al., 2008; Fachri et al., 2011, 2016; Qu and Tveranger, 2016). Our fault facies are based on 3D stochastic deformation band networks that are generated for each deformed fault facies to capture realistic distributions of deformation band sizes and orientations. Fault facies modeling parameters are varied to generate different heterogeneity distributions, and the interplay between the modeling parameters and the reservoir performance is analyzed. The effect of upscaling fault damage zone properties is analyzed by comparing fluid flow simulation results of fine-gridded models and coarse-gridded models. Methodology The investigation of the impact of fault damage zone properties on reservoir-scale fluid flow is based on a conceptual damage zone model of normal faults in porous sandstones, which has been characterized by previous works using extensive outcrop data (Schueller et al., 2013; Qu and Tveranger, 2016). Figure 1 shows a typical deformation-band fault damage zone in porous sandstones and the deformation band arrangement within fault damage zones. The workflow for setting up

The impact of individual deformation bands on fluid flow in sandstones has been extensively studied using physical experiments on core plugs and samples.A review of these studies by Fossen and Bale (2007) shows that, deformation bands, either individually or as clusters, may display a permeability reduction relative to their host rock of up to five to six orders of magnitude.However, the same study concludes that the impact on effective permeability is limited, and controlled by the cumulative thickness of deformation bands and their 3D connectivity.
Understanding the influence of deformation-band fault zones on a reservoir scale requires proper implementa-tion of these features in reservoir models.In general, fault zones are conceived as consisting of a fault core accommodating the bulk displacement and a surrounding damage zone (Caine et al., 1996;Shipton et al., 2005).Conventional reservoir models use transmissibility multipliers to capture the effect of faults zones on cross-fault fluid flow.These are based on the calculated fault zone permeability, fault zone thickness, and permeability of the footwall and hanging-wall grid cells.Fault zone thickness is estimated from correlations against fault throw.The calculation of fault zone permeability is based on an empirical relationship between shale content and displacement (Manzocchi et al., 1999).The "fault zone thickness" in this context is defined as either the separation between the outermost slip surfaces minus the thickness of undeformed lenses, or the thickness of the slip-surface itself in a lacuna (Manzocchi et al., 1999).Identifying the "outmost slip surfaces" is an ambiguous exercise even in well-exposed outcrops.In practice, fault zone thickness is handled as an approximation of the fault core thickness and envisaged as a homogeneous volume of rock representing the deformed rocks present in the fault zone.It is difficult to make a case for also including the damage zone in fault zone thickness because damage zones display clear spatial trends in deformation intensity and commonly are anything but homogenized.Damage zones are up to several orders of magnitude wider than the fault core and extend tens of meters away from the core (Shipton et al., 2002(Shipton et al., , 2005;;Torabi and Berg, 2011).For common reservoir models using 10-100 m XY resolution, damage zone features, if implemented, can be expected to influence fluid flow patterns.Some researchers captured the effect of damage zones by reducing the transmissibility multipliers of the nearby fault connections (Manzocchi et al., 2008).However, the resulting reservoir simulation models do not render explicit fault damage zone heterogeneity and accurate fluid flow calculations in the damage zones.
Several studies have attempted to capture the effect of deformation-band damage zones on fluid flow using reservoir models with explicit representation of damage zones.Some of these studies focus on deriving bulk permeability values for implementation in low-resolution simulation models (Odling et al., 2004;Sternlof et al., 2004;Harris et al., 2007;Fachri et al., 2013b); hence, the spatial scale of their models is generally limited.Similar studies with larger model dimension (Rotevatn et al., 2009;Fachri et al., 2013a) are based on field observation in deriving deterministic upscaled damage zone permeability, leaving room for improvement by using stochastic method where input (e.g., deformation band orientation and volumetric proportion) can be modulated.It has not been clear whether, or under which circumstances, the spatial variation of fault zone properties influences regional flow patterns and reservoir behavior or not.The objective of this paper is to (1) demonstrate an improved workflow in damage zone modeling and (2) investigate the contribution of fault damage zone on reservoir performance and the sensitivity of reservoir performance to fault damage zone heterogeneity, thereby to infer implications regarding fault zone property modeling issues in real-field geologic models.
For this study, we use a multiscale modeling workflow based around the use of fault facies (Tveranger et al., 2005;Braathen et al., 2009) to model a conceptual normal fault in porous sandstone with a damage zone consisting of deformation bands (for the previous applications of fault facies, see Fredman et al., 2008;Fachri et al., 2011Fachri et al., , 2016;;Qu and Tveranger, 2016).Our fault facies are based on 3D stochastic deformation band networks that are generated for each deformed fault facies to capture realistic distributions of deformation band sizes and orientations.Fault facies modeling parameters are varied to generate different heterogeneity distributions, and the interplay between the modeling parameters and the reservoir performance is analyzed.The effect of upscaling fault damage zone properties is analyzed by comparing fluid flow simulation results of fine-gridded models and coarse-gridded models.

Methodology
The investigation of the impact of fault damage zone properties on reservoir-scale fluid flow is based on a conceptual damage zone model of normal faults in porous sandstones, which has been characterized by previous works using extensive outcrop data (Schueller et al., 2013;Qu and Tveranger, 2016).Figure 1 shows a typical deformation-band fault damage zone in porous sandstones and the deformation band arrangement within fault damage zones.The workflow for setting up the models for flow simulations is illustrated in Figure 2.

Model configuration
We use two synthetic reservoir models built using the reservoir modeling software RMS (see Figure 3 for model dimen- sions).Both include a single linear fault dipping 70°to the east (x), with a maximum throw of 100 m.In configuration A, the fault measures 1000 m long and dies out northward and southward inside the model domain, allowing fluids to bypass the fault (Figure 3a).In configuration B, the fault intersects the entire model domain (Figure 3b).Both are mimicking fault configurations found in subsurface reservoirs (i.e., isolated faults [fault configuration A] and linked faults or faults compartmentalizing reservoirs [fault configuration B]).Both reservoir models use a grid resolution of 20 × 20 × 20 m and consist of homogeneous sandstone with a porosity of 0.25 and permeability of 1000 mD.The sandstone reservoir is assumed to sandwiched between impermeable shale, which is not explicitly included in the model.

Fault damage zone model: Gridding
In the conventional reservoir modeling workflow for faults, the alteration of petrophysical properties by the presence of deformation bands in the damage zone is usually ignored, and the cells around the fault surface are assigned permeability values similar to those of the surrounding undeformed host rock.In our workflow, the damage zone properties are treated explicitly based on fault facies approach (Tveranger et al., 2005;Fredman et al., 2008;Braathen et al., 2009;Fachri et al., 2011Fachri et al., , 2016;;Qu and Tveranger, 2016).Two columns of cells on each side of the fault are refined to 1 × 4 × 4 m to accommodate small-scale heterogeneities in fault damage zones.The choice of the area encompassing the fault to be refined is based on the width of the fault damage zone, which is related to the fault throw.A fault damage zone width that is realistic for the fault in our example is defined according to a scaling relationship for deformationband damage zones in siliciclastic rocks reported by Schueller et al. (2013); the estimated half-width of the fault damage zone (i.e., the distance from the fault core to the outer boundary of damage zone) corresponding to maximum fault throw of 100 m is approximately 31 m.The choice of grid resolution for the fault damage zone model is based on the scale of properties to be modeled.
Because the property modeling functions in RMS cannot be performed on grids with local grid refinements (LGRs), the refined region is exported to allow separate handling and property modeling (Figure 2b).Because the geometry of the fault in our example is quite simple, straightforward implementation of LGR can be done in RMS.For Interpretation / November 2017 SP43 more complex cases with curved faults, modelers are recommended to use a "create-fault-grid" function implemented in the fault-modeling tool HAVANA (Norwegian Computing Center), which automates the steps mentioned above (Qu et al., 2015).

Fault damage zone model: Fault facies
The spatial distribution of deformation band density (number of bands crossed along scanlines per meter) in the fault damage zone is modeled in the fault damage zone grid (Figure 2c), based on the use of fault facies (Tveranger et al., 2005;Braathen et al., 2009) and established relationships linking fault facies distributions to fault displacement and distance to the fault core (Qu and Tveranger, 2016).Four fault facies are defined by different categories of deformation band density: high (>20∕m), medium (6-20∕m), low (1-5∕m), and undeformed rock (0∕m); labeled facies H, M, L, and U, respectively, based on deformation band density data collected from 106 scanlines across fault damage zones (Schueller et al., 2013).The population of fault facies in the fault damage zone is performed by the truncated Gaussian simulation (TGS) modeling method implemented in RMS, which has been verified to be able to reproduce the distribution patterns of fault facies in deformationband fault damage zones (Fachri et al., 2013b).
Two important parameters for modeling the fault facies using TGS are the variogram range and facies proportion trend (Fachri et al., 2013b;Qu and Tveranger, 2016), which can be estimated based on empirical data from outcrops.However, because good 3D fault outcrops are rare, it is difficult to obtain good empirical field-data sets that allow quantification of variogram ranges, especially in fault-parallel direction; thus, it is a main objective of this paper to investigate the sensitivity of reservoir-scale fluid flow to these modeling parameters.
Five fault facies models for the sensitivity study are generated by modulating the modeling parameters (Figure 4).Fault facies models 1−4 are modeled with the  same facies probability trends but different variogram ranges.Fault facies model 5 is a model for comparison purpose, which uses the same variogram range and the same volumetric fractions of fault facies as scenario 3, but shows random facies distributions without the constraint of facies probability trends.

Fault damage zone model: Petrophysical properties
The fault facies model described above explicitly renders the heterogeneity caused by the variation of deformation band density in fault damage zones at the chosen level of grid resolution.However, before conducting flow simulation on models with different facies distributions, we need to estimate effective permeability for each deformed fault facies.Previous studies calculate effective permeability of fault facies with different deformation band densities by using analytical methods (arithmetic or harmonic averages), assuming that the bands are parallel to each other (Rotevatn et al., 2009;Qu and Tveranger, 2016).These estimates do not apply to the conjugate patterns of deformation bands (Figure 1b; Shipton andCowie, 2001, 2003;Fossen et al., 2007;Brandes and Tanner, 2012;Imber et al., 2012).In this paper, we present an alternative method for estimating the effective permeability of deformed fault facies: flow-based upscaling of 3D models in which the conjugate arrangement of deformation bands can be honored.This is presented in detail in the "Deformation band models" section.

Merging fault damage zone model into original reservoir model
After the refined grid is populated with petrophysical values of deformed rocks, it is merged with the original global model (Figure 2e).The merging step is performed using HAVANA.Flow simulations are performed on the merged models with detailed fault damage zone permeability structures to study the influence of fault damage zone properties.
There is an option to upscale the high-resolution fault damage zone properties models (Figure 2d) to a coarser fault-zone grid (Figure 2f) before merging into the original reservoir model.This reduces the number of cells in the simulation model, but also decreases the level of geologic detail.The effect of upscaling is also investigated in this paper.

Deformation band models
Except for facies U, which is composed of undeformed host rock, the other three facies (H, M, and L) include networks of deformation bands in a host-rock matrix.To obtain the effective permeabilities for these, we generate a series of high-resolution models that explicitly render the deformation-band networks and rock matrix.The models are conceptual, with the modeling parameters based on published outcrop studies of deformation bands in fault damage zones (Fossen and Hesthammer, 1997;Shipton andCowie, 2001, 2003; Fos-    , 2007;Brandes and Tanner, 2012;Imber et al., 2012).The process involves five steps and is illustrated in Figure 5.
Step 1 is stochastic modeling of deformation band networks, performed by a tool for stochastic modeling of subseismic-scale faults (HAVANA; Hollund et al., 2002).The input modeling parameters are listed in Table 1, which includes orientation and size distributions of deformation bands, length-height ratios of deformation bands, simulation volume, and the number of bands to be simulated.The strike and dip set up for deformation bands are based on outcrop observations.In plan view, deformation bands are not strictly parallel to each other, which defines thin, lozenge-shaped blocks of relatively undeformed sediment between deformation bands (Figure 1a; Shipton andCowie, 2001, 2003).In cross sections, deformation bands within the damage zones dip synthetically and antithetically to the main fault, forming a conjugate geometry (Figure 1b; Shipton andCowie, 2001, 2003;Fossen and Bale, 2007;Fossen et al., 2007;Brandes and Tanner, 2012;Imber et al., 2012).The output is the surface point data, in which each surface represents a band.Multiple realizations of deformation band networks are generated for each type of fault facies (H, M, and L).The generated surface point data are treated as input fault surface data for building structural models in RMS (step 2).
The structural models are then used for gridding (step 3).A volume of 1 × 3 × 3 m within the structural model is gridded into one million cells with a resolution of 1 × 3 × 3 cm, with the deformation bands (or "faults") stair stepped (Figure 6).The choice of the volume size to be gridded is a compromise between the two: a higher grid resolution to honor more realistic band geometries and less total number of cells that the upscaling program can handle within a reasonable time.
Deformation bands are defined by cells adjacent to stair-stepped faults, and permeability values are assigned to the deformation bands and matrix (step 4, Figure 6).Note that there are always two columns of cells adjacent to each fault, and only one column of faulted cells is defined for each deformation band.This is done by using internal programming language of RMS.Thus, the deformation band in the model is one cell thick (1 cm).
Step 5 is to obtain the upscaled permeability of the 1 × 3 × 3 m volume.Single-phase flow-based upscaling was carried out in RMS.A sealed boundary condition, which has been used for damage zone permeability upscaling by previous authors (Odling et al., 2004;Harris et al., 2007), was applied here.The input is the fine-scale grid model from step 4, in which deformation bands and matrix are explicitly represented and assigned representative permeability.For simplicity, the permeability of the deformation bands (k db ) and matrix (k matrix ) was assumed to be uniform and isotropic: 1000 mD for the matrix and three to five orders magnitude lower (1, 0.1, and 0.01 mD) for the deformation bands.For each model, the outputs are upscaled permeability values in the three directions (x, y, and z).Table 2 lists the minimum, maximum, and average values of the upscaled permeabilities of multiple realizations of fault facies.Coefficient of variation (Cv) of the upscaled permeability values of multiple realizations for each fault facies is also calculated to look at the representativeness of the upscaled permeability values (Corbett and Jensen, 1992;Nordahl et al., 2005).It shows the variation of the upscaled permeability values of multiple realizations of facies H and L is quite small (Cv < 0.5), whereas K Y values (upscaled permeability in fault-parallel direction) of multiple realization of facies M have a larger variation (Cv > 0.5), especially when the permeability contrast between deformation band and undeformed rock is high.Thus, the upscaled permeability of facies M should be used with caution for permeability population in grid models with cells of much larger scale.In this paper, these permeability values of the 1 × 3 × 3 m volume are populated in grid cells of similar size.

Fluid flow simulation models
Fluid flow simulations were carried out to study the effect of fault damage zone and the effect of changing specific input parameters used in the fault damage zone model (variogram range and facies probability trends) on reservoir performance.Models run for flow simulations are listed in Table 3.The flow impact of the fault damage zone when combined with different fault core properties was investigated by modulating the fault core transmissibility multiplier (FTM) values.The effect of upscaling fault damage zone properties is analyzed by comparing fluid flow simulation results of models prior to and after upscaling.
Grid resolution of the simulation models prior to upscaling is 1 × 4 × 4 m in the fault damage zone region, 20 × 20 × 20 m in the rest of the model (Figure 2d and  2e).Grid resolution of the simulation models after upscaling is 20 × 20 × 20 m in the entire model (Figure 2f and 2g).The grid type is a corner point grid.The total numbers of cells in the simulation models are listed in Table 3.
Based on the method for calculating effective permeability presented in the previous section, each fault facies in the fault damage zone model was assigned the mean value of their respective permeability distribution (Table 2).Previous calculations show that although the porosity of deformation bands may be more than one order of magnitude lower than the surrounding host rock, this does not significantly alter the bulk porosity of the rock containing deformation bands.There is very little difference between the bulk porosity of rocks with different deformation band density values (Antonellini and Aydin, 1994;Qu and Tveranger, 2016).In this work, a porosity of 0.25 was assigned to the undeformed rock, and 0.23 to the deformed fault facies.
One vertical water injector was placed in the hanging wall (360 m from the fault), and one vertical producer was placed in the footwall (360 m from the fault).Fluid flow simulations were conducted using Eclipse.There was no flow across model boundaries.All layers were perforated.Wells were controlled by predefined rate targets with pressure limits.All simulations were run until the water cut in the production well exceeded 90%, but with a maximum limit of 20 years.Details of rock, fluid, and production parameters are listed in Table 4, in which the material constants are derived from published literature (Fachri et al., 2011).Initial saturation and pressure distribution are defined by the simulator using a user-defined depth for the oil-water contact (which is in our case set to below the bottom of the reservoir), fluid pressure at a reference depth, and fluid pressure gradient.Uniform relative permeability curves were assigned to all the cells in the model (Qu and Tveranger, 2016).

Results and discussion
The effect of fault damage zones in combination with different fault core transmissibility multipliers  The effect of the fault damage zone on fluid flow is decreased when a baffling fault core (i.e., a lower FTM value) is present.Figure 7b and 7c shows that the difference in oil recovery and time to water breakthrough between the damage zone present models and damage zone absent models is decreasing with the decrease of FTM value from one to zero.The damage zone-present models and the damage zone-absent models even show the same reservoir responses when an FTM of zero is applied (i.e., when the fault core is completely sealing; Figure 7b and 7c).Solely based on this, one might easily conclude that the presence of a fault damage zone can be ignored in a noncompartmentalized reservoir, if the fault core permeability is sufficiently low.However, our results suggest that the presence of a fault damage zone appears to affect flow patterns in these instances because unswept oil can be observed in the footwall damage zone (Figure 9).This effect cannot be modeled by the implicit modeling method, which uses a reduced fault transmissibility multiplier value to capture the flow effect of a baffling fault damage zone; as Figure 9f shows, the effect of unswept oil in the down-flow side of the fault is not captured by the FTM-only model even though an extremely low value of zero is used.If the   fault damage zone volume is large, the oil volumes retained in the damage zone could be substantial and require an adaption of production strategies to ensure their recovery.This would necessitate the use of a model with an explicitly rendered damage zone.
The results also show that reservoir performance is not sensitive to the variation of an FTM in damage zonepresent models (Figure 7).

Configuration B: Fault intersects entire model domain
In configuration B, replicating a compartmentalized reservoir in which flow is forced to move through the fault, the inclusion of a fault damage zone generally has a negative effect on production.Damage zone-present models show lower oil recovery and delayed water breakthrough compared with the damage zone-absent models (Figure 7d-7f).Cross-fault flow is inhibited, and less area is swept due to the introduction of a damage zone (Figures 10 and 11).
The influence of a fault damage zone on reservoir performance is significant even in the presence of a low-transmissibility fault core (FTM ¼ 0.01, Figure 7d  and 7e).When an FTM of 0.01 is combined, oil recovery after 600 days was up to 0.5 in the damage zone-absent model, whereas it was less than 0.003 in the damage zone-present models (Figure 7d).Water breakthrough happened after 500 days of production simulation in the damage zone-absent model, whereas it did not happen during the maximum simulation time of 20 years in the damage zone-present models (Figure 7f).No flow simulation using an FTM of zero was performed for configuration B, because this setting would effectively isolate the footwall and hanging wall from each other.
In total, not including fault damage zone properties in configuration B reservoir models, even in the presence of a low-transmissibility fault core (FTM ¼ 0.01), can lead to erroneous flow simulation results.

The effect of fault facies modeling parameters
In configuration A, reservoir performance is not sensitive to the fault facies modeling parameters.Models with different variogram ranges and facies probability trends show similar reservoir responses (Figures 7a-7c, 12a, and 12b) and flow patterns (Figures 8 and 9).
In configuration B, reservoir performance is affected by varying fault facies modeling parameters.The plots of reservoir response parameters (water cut and oil recovery after 4000 days) versus models with different variogram ranges show that the reservoir response is more sensitive to the variogram range in the faultperpendicular direction than in the fault-strike and faultdip directions (Figure 12e and 12f).Models with a larger variogram range in the fault-perpendicular direction (10 m in the test models) show much higher oil recovery and higher water cut than the models with a smaller variogram range (2 m) in the fault-perpendicular direction (Figure 12e and 12f).Regional flow patterns are also affected by the variogram range (Figure 10a-10d).An irregular water front in homogeneous host rock is observed (Figure 10d).Once the fluid finds a way to cross the fault zone, it quickly spreads out into the host rock, leaving some part of the fault zone unswept (Figure 10c).The sensitivity of the reservoir response to the variogram range is decreased in the presence of a low FTM (FTM ¼ 0.01); all four models with different variogram ranges show very low oil recovery (Figure 12f), and water is hardly produced (Figure 12e).
The reservoir responses of configuration B models are quite sensitive to facies probability trends, even in the presence of a low FTM. Figure 12c and 12d shows the plots of oil recovery and water cut versus models with and without the constraint of facies probability trends (i.e., fault facies models 3 and 5 in Figure 4).The model without the constraint of facies probability trends shows much higher oil recovery and much higher water cut than the model with the constraint of facies probability trends (Figure 12c and 12d).In the presence of a low FTM, the reservoir model with random fault facies distributions has been largely swept at the end of the simulation (Figure 11e), compared with the little swept regions in the models constrained by facies probability trends (Figure 11a-11d).Note that the two models for comparison are modeled with the same variogram range, and the volumetric fractions of fault facies in these two produced models are the same.The only reason for the difference in reservoir performance is the distribution of fault facies.This highlights the significance of fault damage zone heterogeneity to reservoir performance and the importance of characterizing and modeling fault facies distributions.
These sensitivity tests show that the fine-scale permeability structure of fault damage zones has a different degree of significance to configuration A models and configuration B models.This finding gives a hint to whether or not upscaling the fault damage zone would affect reservoir performance by losing permeability structures, which is addressed in the "Influence of upscaling" section below.

Influence of upscaling
Results show that the reservoir performance of configuration A models does not seem to be affected by the upscaling (Figure 13a and 13b).The models prior to and after upscaling display the same oil recovery (Figure 13a) and a similar water cut (Figure 13b).This is consistent with the finding above, which shows that the reservoir performance of configuration A models is not sensitive to fault facies distributions and there is no loss to average up the permeability structures inside the damage zone.
In configuration B, contrasting reservoir responses between the fine-scale models and the upscaled models are observed (Figure 13c and 13d).The oil recovery and water cut are reduced after upscaling due to the smoothing up of some permeability structures critical to fluid flow.This is also consistent with the finding above that the reservoir performance of configuration B models are quite sensitive to fault facies distribtuions.In this case, using a high-resolution fault zone grid in simulation models seems a necessity, but it comes with a significant CPU cost.In the conducted test, the models were scaled up 20 times in the fault-perpendicular direction, and sealed boundary conditions were applied.It will be future work to check if a limited upscaling and a different boundary condition could yield a reliable flow simulation result.

Implications on handling of fault damage zones in reservoir models
The overall impact including a fault damage zone has on reservoir performance can be attributed to the fact that the fault damage zone is a 3D volume, which can accommodate and retain oil and injected water, and therefore it influences waterfront movement in the region outside the fault damage zone.This study highlights the potential importance of including explicitly rendered the damage zones in reservoir models.
This paper finds that the noncompartmentalized reservoirs in which flow can move around the fault tips have a lower requirement for precise fine-scale fault damage zone property modeling and the upscaling procedure does not affect the reservoir performance.In this case, the handling of fault damage zones in simulation models can be easier: Upscaled properties can be populated in a coarser grid directly.For example, the cell-permeability values in the upscaled damage zone model shown in Figure 2f can be populated elsewhere of a similar case.This skips the fine-scale damage zone property modeling, but a thorough database of effective petrophysical properties of different geologic scenarios is a prerequisite.Effort going in this direction can be found in Fachri et al. (2013b).
The compartmentalized reservoirs require much more precise fine-scale fault damage zone property modeling.A prerequisite to reach this goal is a systematic data set from outcrops, which can provide input for fault facies modeling.Apart from that, seismic characterization can be a supplement to condition the modeling, and forward seismic studies have shown promising results toward this direction (see Botter et al., 2017;Kolyukhin et al., 2017).

Conclusion
We have investigated the impact of incorporating a field-scale 3D fault damage zone on simulated fluid flow using two scenarios: one in which fluid is allowed to flow around the fault zone (configuration A) and one in which fluid is forced to move through it (configuration B).The fault facies model of the fault damage zone is based on empirical data collected from outcrops.Flow-simulation testing was carried out using models, in which the facies variogram ranges, facies probability trends, grid resolution, and FTM were varied.Simulation results from all models were compared to establish the impact of the different fault damage zone model input parameters on simulated reservoir behavior.
In configuration A, in which fluid is allowed to flow around the fault tips, the deformation-band fault damage zone generally shows a positive effect to production; this effect decreases when the FTM decreases.The fault damage zone will contain undrained oil when the fault core is modeled as completely sealing.This implies the necessity of representing the fault damage zone explicitly in reservoir models, especially if the fault damage zone volume is suspected to be large.Reservoir performance in configuration A is not sensitive to fault damage zone heterogeneity.The upscaled models, in which the fault damage zone is represented by one to two columns of cells, produce the same reservoir response as the initial high-resolution models.
In configuration B, in which fluid is forced to cross the fault, the presence of a fault damage zone has a neg-ative effect to production, and this effect is still significant in the presence of a low-transmissibility fault core.Reservoir performance is generally influenced by varying the fault damage zone heterogeneity, and using upscaled fault damage zone properties at the scale used here (20 × 20 × 20 m) did not reproduce reservoir behavior of the initial fine-scale model.However, a less radical upscaling may achieve a better balance between producing reliable forecasts of reservoir behavior and the required number of cells in simulation models.
This study highlights the importance of including fault damage zones in subsurface reservoir models.This calls not only for improved tools and workflows, but most importantly, it also calls for new empirical data set that will allow improved 3D characterization of fault zones for reservoir modeling purposes.

Figure 1 .
Figure 1.Illustration of typical deformation-band fault damage zones in porous sandstones.(a) Fault damage zone of Bartlett Fault, Utah.(b) Schematic illustration of the typical conjugate arrangement of deformation bands in the damage zone around a simple fault, modified from Fossen and Bale (2007).

Figure 2 .
Figure 2. Illustration of workflow performed in this paper.(a) Synthetic reservoir model containing an east-dipping fault.The maximum fault throw is 100 m.(b) Fault-zone grid.(c) Fault facies model and deformation band network models for each fault facies.(d) Fault damage zone permeability model.(e) Merged model.(f) Upscaled fault damage zone permeability model.(g) Merged model with upscaled fault damage zone permeability.

Figure 3 .
Figure 3. Two model configurations used in this paper.(a) Configuration A: The fault dies out inside the mode domain.(b) Configuration B: The fault intersects the entire model domain.

Figure 5 .
Figure 5. Flowchart for modeling deformation band network and obtaining the upscaled permeability of fault facies.
the mean direction parallel to the fault and standard deviation of 10°D ip Set 1: synthetic to the main fault, mean dip angle = 70°, standard deviation = 5°S et 2: antithetic to the main fault, mean dip angle = 70°, standard deviation = 5°L ength (L) D max ¼ 0.0016L 0.5 (Fossen and Hesthammer, 1997) Displacement range (D): 0.003-0.1 m Length-to-height aspect ratio 2 Number of simulated bands 110, 75, and 40 for facies with high, medium, and low deformation band densities Dimension of simulation volume 3 × 10 × 7 m Interpretation / November 2017 SP45 sen et al.

Figure 6 .
Figure 6.Deformation band network models for (a) facies H, (b) facies M, and (c) facies L. The left panel shows deformation band network models, where deformation bands are represented as fault surfaces in structural models.The structural model covers a range of 3 × 7 × 10 m.A volume of 1 × 3 × 3 m from the center of the structural model is gridded into one million cells.The middle panel shows models in which deformation bands (green) and matrix (red) are explicitly represented.The right panel shows cells intersected by deformation bands.
Configuration A: Fault dies out inside model domainIn configuration A, in which fluid is allowed to flow around the fault tips, the presence of a deformationband fault damage zone generally has a positive effect on production because it forces fluids to flow around the fault tips, which improves sweep and delays the time to water breakthrough (Figures7a-7c and 8).

Figure 7 .
Figure 7.Comparison of reservoir responses after a certain simulation time versus models with and without damage zones.The displayed damage zones present models in each plot are modeled with the same facies probability trends and different variogram ranges.(a) Oil recovery after 1000 days of configuration A modes with and without damage zones.FTM ¼ 1. (b and c) Oil recovery after 1000 days and time to water breakthrough of configuration A models with and without damage zones in combination with different FTM values.(d) Oil recovery after 600 days of configuration B models with and without damage zones.FTM ¼ 1. (e and f) Oil recover after 600 days and time to water breakthrough of configuration B models with and without damage zones in combination with different FTM values.(f) Note that water breakthrough did not happen in configuration B models in the presence of the damage zone and FTM of 0.01 during the maximum simulation time of 20 years.

Figure 8 .
Figure 8.Oil saturation of configuration A models after two years of production simulation.FTM ¼ 1. (a-d) Damage zone-present models with the constraint of the same facies probability trends and different variogram ranges.(e) Damage zone-present model without the constraint of facies probability trends.k matrix ∕k db ratio is 10 4 .(f) Model with only fault transmissibility multiplier applied and the damage zone is omitted.Layer 5 and crossing-well section are displayed in each panel.

Figure 9 .
Figure 9. Oil saturation of configuration A models at the end of production simulation.FTM ¼ 0. (a-d) Damage zone-present models with the constraint of the same facies probability trends and different variogram ranges.(e) Damage zone-present model without the constraint of facies probability trends.k matrix ∕k db ratio is 10 4 .(f) Model with only fault transmissibility multiplier applied and the damage zone is omitted.Footwall block and crossing-well section are displayed in each panel.

Figure 10 .
Figure 10.Screenshots of oil saturation of configuration B models after 10 years of production simulation.FTM ¼ 1. k matrix ∕k db ratio is 10 4 .The middle layer (layer 5) and the crossing-well section are displayed in each panel.(a-d) Damage zone-present models with the constraint of the same facies probability trends and different variogram ranges.(e) Damage zone-present model without the constraint of facies probability trends.(f) Model without damage zone.Displayed time for the model without damage zone is 15 months, when production is completed.

Figure 11 .
Figure 11.Oil saturation of configuration B models in the end of production simulation.FTM ¼ 0.01.(a-d) Damage zone-present models with the constraint of the same facies probability trends and different variogram ranges.(e) Damage zone-present model without the constraint of facies probability trends.(f) Model with fault damage zone omitted.Fault transmissibility multiplier of 0.01 is applied to all the models.Layer 5 and crossing-well section are displayed in each panel.k matrix ∕k db ratio is 10 4 .

Figure 12 .
Figure 12.Comparison of reservoir responses of models with different fault facies modeling parameters.(a and b) Oil recovery and water cut of configuration A models with and without the constraint of facies probability trends.(c and d) Oil recovery and water cut of configuration B model with and without the constraint of facies probability trends.(e and f) Oil recovery and water cut of configuration B models with different variogram ranges.

Figure 13 .
Figure 13.Comparison of reservoir responses of models prior to and after upscaling.(a and b) Oil recovery and water cut of configuration A models prior to and after upscaling.(c and d) Oil recovery and water cut of configuration B models prior to and after upscaling.

Table 1 .
Input parameters for modeling deformation band network.

Table 2 .
Upscaled permeabilities of deformation band network models for different fault facies (k db = permeability of deformation band, k matrix = permeability of undeformed rock, K = upscaled permeability, and Cv = coefficient of variation).

Table 4 .
Dynamic properties used for fluid flow simulations.