## Abstract

It is well established that compaction bands (CBs), joints and faults are often present in the same rock volume in the Jurassic aeolian Aztec Sandstone, exposed in the Valley of Fire State Park, Nevada, USA. Because the permeability of CBs can be one or more orders of magnitude less than the matrix permeability, and joint permeability, depending on its aperture, can be several orders of magnitude greater than matrix permeability, the combined effect of these structures on subsurface flow can be complex and substantial. In this study, we investigate the effects of a variety of intersecting geological structures on fluid flow. This is accomplished by performing two- (2D) and three-dimensional (3D) permeability upscaling and waterflood simulations over areas/volumes populated by hydraulically interacting geological features. The regions considered are approximately the size of typical grid blocks used for reservoir or aquifer flow simulations, so the results are of practical interest. The systems studied include models with two sets of vertical CBs intersecting at various angles, an inclined CB set intersecting a vertical CB set, a joint set intersecting a CB set at various angles, and a small fault and its damage zone overprinting a CB set. Our numerical results quantify the impact of these composite structures on subsurface flow and show, for example, that the intersection angle of two sets of structures can have a considerable effect on the upscaled directional permeability. In addition, waterflood simulations demonstrate that the efficiency of oil recovery can be significantly impacted by the direction of flow relative to the orientation of intersecting geological structures.

Structural heterogeneities, such as dilatant fractures (joints), shear fractures (faults) and contraction/compaction structures (compaction and shear bands and solution seams), are widely observed in clastic rocks. Each of these structure types occur in a specific geological and geomechanical environment and have a specific genesis. Each type has its own geometry, spacing, distribution, connectivity and hydraulic properties, which either enhance or impede subsurface fluid flow. Permeabilities of these structures may, on average, be a few orders of magnitude higher or lower than those of the corresponding matrix rocks. Therefore, structural heterogeneities are essential elements of subsurface fluid flow, and should be included in large-scale flow simulations of aquifers, reservoirs and basins (Aydin 2000).

In recent years, researchers have extensively studied the influence of shear bands, compaction bands (CBs), joints and faults on subsurface fluid flow, although generally only one feature at a time was considered (e.g. Pittman 1981; Antonellini & Aydin 1994; Fowles & Burley 1994; Matthai *et al.* 1998; Taylor *et al.* 1999; Flodin *et al.* 2001; Main *et al.* 2000, 2001; Jourde *et al.* 2002; Shipton *et al.* 2002; Myers & Aydin 2004; Sternlof *et al.* 2004, 2006; Ahmadov *et al.* 2007; Fossen & Bale 2007; Kolyukhin *et al.* 2010). Matthai *et al.* (1998) and Jourde *et al.* (2002) noted that faults with highly permeable slip planes are the most dominant component because they are more continuous and thus influence flow at the reservoir scale. High-permeability joints with some degree of connection were found to be the next most important element, while deformation bands (both shear and compaction types) with somewhat large spacing had only a relatively minor impact on flow.

Taylor *et al.* (1999) combined field mapping and numerical modelling of palaeo-flow systems in sedimentary rocks and found that the distribution of flow within conductive joint sets is influenced by the geometric pattern of joints and the hydraulic properties of both joints and matrix. Using a reservoir characterization approach and extensive core coverage, Shipton *et al.* (2002) identified that faults in aeolian sandstone consist of three types of components: deformation bands (primarily shear bands); the fault core; and slip surfaces. These authors performed sensitivity analyses, which indicated that the effective transverse fault-zone permeability is most dependent on the thickness of the fault core.

Using a geologically-based methodology, Flodin *et al.* (2001) demonstrated the impact of variable fault-zone permeability on reservoir-scale flow models. Sternlof *et al.* (2004) numerically studied the influence of compaction band (CB) pattern on effective sandstone permeability at small scales (*c.* 1–10 m) and demonstrated the potential for compaction band arrays to induce one–three order-of-magnitude effective permeability reduction in sandstone. Sternlof *et al.* (2006) modelled flow at larger scales using mapped CB distributions, and noted the potential significance of CBs in aquifer- and reservoir-scale simulations. Fossen & Bale (2007) presented mathematical calculations indicating that very high permeability contrasts between deformation bands and rock matrix and/or exceptionally high band concentrations are required for deformation bands to significantly affect production rates. Using an upscaling procedure, Kolyukhin *et al.* (2010) demonstrated that the effect of deformation bands (primarily shear bands) on permeability in the fault damage zone depends on band clustering, orientation and the distribution of permeability anomalies or ‘holes’ in the bands.

With the exception of Matthai *et al.* (1998) and, to a limited extent, Taylor *et al.* (1999) and Shipton *et al.* (2002), most of the previous work only considered the influence of a single type of geological structure on subsurface fluid flow. This is because they focused on the formation, mechanics and flow properties of a particular structure type. Consequently, they isolated the structure of interest in order to analyse its behaviour. However, in many locations worldwide, and in the Valley of Fire State Park (NV) (see Fig. 1a inset for location) in particular, different types of geological structures are commonly present together in the same rock volume. CBs and shear bands, joints and faults (Fig. 1a–d) overprint each other within the Aztec Sandstone exposed throughout the park, which has been the focus of several doctoral dissertations (Antonellini 1994; Myers 1999; Taylor 1999; Flodin 2004; Sternlof 2006) and postdoctoral projects (Jourde *et al.* 2002; de Joussineau & Aydin 2007; de Joussineau *et al.* 2007) at Stanford. For example, Figure 1b displays two sets of CBs, one at a high angle to bedding and the other subparallel, and Figure 1c shows a set of joints intersecting a high-angle set of CBs with an intersection angle of about 30°. In Figure 1d, there is a small fault (underneath the measuring tape in the picture) which is subparallel to the CBs and overlaps with one CB. In addition, a series of asymmetric splay joints in the damage zone of the small fault intersect the CB set at a low angle of about 30°.

As the structures described above have contrasting permeabilities, their occurrence within the same rock volume will, in general, result in very different effects on subsurface flow compared to the effect of a single type of geological structure. It is therefore evident that overlapping patterns of different types of structures deserve more attention. In this work we focus on the patterns observed in the Valley of Fire State Park, Nevada, which is located about 60 km NE of the city of Las Vegas. The predominant outcrop unit in the park is the Jurassic Aztec Sandstone, with spectacular colour patterns that were diagenetically controlled (Eichhubl *et al.* 2004). The Aztec Sandstone has a stratigraphic thickness of about 1400 m in the park and its vicinity (Bohannon 1983; Marzolf 1990). Texturally, it is a well-sorted quartz arenite with a mean grain size of about 0.25 mm, a porosity range of 15–25% and a permeability range of 100–5900 mD (millidarcies, where 1 mD=10^{−15} m^{2}) (Antonellini & Aydin 1994; Antonellini *et al.* 1994; Flodin *et al.* 2004).

In this work we use numerical models to investigate the fluid-flow properties of the following structure configurations: (1) two sets of vertical CBs with various intersection angles; (2) an inclined CB set with various inclination angles intersecting a vertical CB set; (3) joints overprinting and intersecting previously existing CBs at various angles; and (4) a small fault with an asymmetric damage zone overprinting and intersecting earlier CBs. We present numerical results for equivalent (upscaled or effective) permeability and for waterflood response to demonstrate the impact of intersecting geological structures on fluid flow. The regions considered are about the size of grid blocks used in large-scale reservoir or aquifer flow simulations, so our findings are practically relevant. We consider a range of structural parameters (e.g. intersection angles) and various CB permeability values. The results clearly demonstrate the importance of treating combined configurations rather than individual patterns in isolation.

## Fluid-flow simulation model

Both single-phase and two-phase (oil–water) flow simulations are performed in this work. The single-phase flow equation is given by:
(1)
where *p* is fluid pressure, *Q* represents sources and/or sinks (e.g. injection/production wells), μ is fluid viscosity and *k* is the absolute rock permeability. In this paper, *k* is assumed to be locally isotropic in all features (e.g. rock matrix, CBs, joints).

Two-phase (oil–water) flow is governed by the following equations (gravitational and compressibility effects are neglected for convenience):
(2)
(3)
where ϕ is rock porosity, the subscripts w and o indicate the water phase and oil phase, respectively, and *S* and *k*_{r} denote the saturation and relative permeability for each phase. The relative permeabilities, *k*_{rw} and *k*_{ro}, are functions of *S* and are typically determined from core experiments. In this work we neglect capillary pressure so *p*_{w}=*p*_{o}. The system is closed by requiring *S*_{w}+*S*_{o}=1.

As described in Sternlof *et al.* (2006) and Ahmadov *et al.* (2007), a specialized finite-volume discretization technique, developed by Karimi-Fard *et al.* (2004), can be applied in the flow simulations. Using this approach, thin features such as fractures and compaction bands are represented as linear elements in unstructured 2D models and as planar elements in 3D models. This simplifies grid generation and reduces the total number of cells in the model. The cell-to-cell transmissibilities obtained from this discretization technique, along with the cell geometrical information, are input to Stanford's General Purpose Research Simulator (GPRS) for fluid-flow simulation. Detailed information on GPRS can be found in Cao (2002). Note that, in the simulations performed here, we use a two-point flux approximation, so some inaccuracy can result if the grid is significantly (and systematically) nonorthogonal. This is not expected to be of significant concern for the types of models considered in this study.

We quantify the effects of the various geological structures on fluid flow by computing the upscaled (or equivalent) permeability over the model domain. The model domains considered are typically of about the size of a single grid block (e.g. 20 m in the *x*- and *y*-dimensions) in a large-scale flow model. Thus, the upscaled values reported correspond to grid-block permeabilities for blocks containing various configurations of features.

We describe our permeability upscaling procedure with reference to Figure 2, which depicts a 2D rectangular region of physical dimensions *L*_{x} and *L*_{y}. In this region, permeability varies from cell to cell. We set the viscosity μ =1 (this does not impact the results for upscaled permeability) and proceed by solving two separate single-phase flow problems, with imposed pressure difference (Δ*p*) in the *x*- and *y*-directions, respectively. No-flow conditions are specified on the other boundaries. For each problem, we compute the total flow rate through the system *Q*. Designating the prescribed pressure differences as Δ*p*_{x} and Δ*p*_{y}, and the corresponding flow rates as *Q*_{x} and *Q*_{y}, upscaled permeability components in the *x*- and *y*-directions ( and ) can be calculated directly using the outlet averaging method (Jourde *et al.* 2002; Flodin *et al.* 2004):
(4)

These upscaled permeability values suffice for many applications, but they define permeability components only in the co-ordinate directions and not a proper full-tensor quantity. When full-tensor effects are important, as can be the case in systems with oriented features, an alternate treatment must be applied. Here we use the area-averaging approach of Wen *et al.* (2003). In this procedure, the average Darcy velocity 〈**u**〉 is related to the average pressure gradient 〈∇*p*〉 in terms of the upscaled permeability tensor **k*** through the expression 〈**u**〉=−**k*** · 〈∇*p*〉. The quantities 〈**u**〉 and 〈∇*p*〉 are computed by averaging over the entire domain as follows:
(5)
(6)
where *j*=1 designates the problem with the prescribed pressure drop in the *x*-direction and *j*=2 the problem with the pressure drop in the *y*-direction, and *A*=*L*_{x}*L*_{y} is the area of the target domain. After obtaining and (*j*=1, 2), the following linear system is solved (using a least-squares method as the system is overdetermined) to provide the components of the full-tensor, **k***:
(7)

Note that the last equation is used to enforce symmetry in **k***. See Wen *et al.* (2003) for full details and Li *et al.* (2011) for an alternate approach for maintaining symmetry in **k***.

Analogous treatments can be used for 3D problems. In this case, we need to solve three flow problems with pressure drop in the *x*-, *y*- and *z*-directions. In the 3D results reported here, rather than compute the six components required to characterize a full-tensor, **k***, we apply only the outlet averaging method (Equation 4). This provides the three permeability components parallel to the co-ordinate directions.

In addition to computing upscaled permeability, we also perform waterflood simulations for several cases. This enables us to examine the effects of various types of geological structures on oil recovery. For these simulations we solve Equations (2) and (3), subject to various boundary conditions, using GPRS.

## Parameters used in fluid-flow simulations

The parameters used for our flow simulations (Table 1) are based on previous work on the Aztec Sandstone (Flodin *et al.* 2001; Jourde *et al.* 2002; Sternlof 2006; Ahmadov *et al.* 2007). Through Lattice–Boltzmann flow simulations, Sternlof (2006) concluded that permeability reduction due to compaction within the bands is at least two orders of magnitude. Aydin & Ahmadov (2009) and Sun *et al.* (2011) also used computational methods and reported that the permeability of CBs is about one order of magnitude lower than that of the rock matrix. Using an image-processing method, Torabi & Fossen (2009) concluded that the permeability range of deformation bands can be a factor of one–four orders of magnitude less than that of the surrounding matrix.

Based on these reported values, for our base-case computations, we assign the permeability of CBs to be two orders of magnitude less than that of the rock matrix. This corresponds approximately to the middle of the range of permeabilities reported for deformation bands in general. We also use other values for CB permeability (specifically 10^{−1} and 10^{−3} of the matrix permeability) in some of our computations to study the impact of this parameter on fluid flow. The permeability of both joints and slip surfaces is calculated using a parallel plate model with an equivalent porous media representation (Taylor *et al.* 1999). We specify an aperture of 0.1 mm for joints and slip surfaces. Using the parallel plate model, this gives a corresponding permeability of about 10^{6} mD (in this work we specify a permeability of 10^{6} mD for joints and slip surfaces).

To simulate two-phase (oil–water) flow, additional data are required. Table 2 provides the relative permeabilities for oil and water. This data, used by Sternlof *et al.* (2006), is representative of a light crude oil (non-wetting phase) and formation water (wetting phase) in a granular medium similar to the Aztec Standstone. The densities of the oil and water are set to 49 and 63 lbm/ft^{3}, respectively, while the viscosities of the oil and water are 2 and 1 cP (centipoise), respectively. In all two-phase flow simulations, it is assumed that the entire rock volume is initially saturated with oil. As stated earlier, capillary pressure is neglected in these simulations.

## Fluid-flow modelling results

In the following, we use the methodology and parameters introduced in the previous sections to investigate the impacts of various types of overprinting and intersecting geological structures on fluid flow in the Aztec Sandstone. Except where otherwise stated, the upscaled permeabilities reported below are computed using the outlet average procedure (Equation 4).

### Two sets of vertical compaction bands

In the Aztec Sandstone, multiple sets of vertical compaction bands (CBs) with different strike directions are commonly present in the same rock volume. We take a set of CBs in an area of 20×20 m from a field map by Sternlof (2006). In our model, we fix the CB set oriented in the *y*-direction and then overprint the other (identical) set on the first set with various intersection angles α. Because the two CB sets would be coincident when α=0°, we let α vary from 5° (subparallel CB sets) to 90° (perpendicular CB sets) in our computations. Figure 3 shows three model configurations with α=5° (upper), 45° (middle) and 90° (lower). These models are discretized for flow simulation using from about 10 000 to 30 000 cells (different numbers of cells are required for the various models because the CB geometries and intersections, which must be resolved by the mesh, differ substantially between models). Rock matrix cells are represented using triangular elements, while CB cells are represented by linear segments.

Figure 4 shows the upscaled permeabilities in the *x*- and *y*-directions, and , normalized by the matrix permeability, *k*_{m}, v. the intersection angle, α, of the two CB sets. The input value of *k*_{CB}/*k*_{m} is taken to be 10^{−1}, 10^{−2} and 10^{−3}, as indicated in each subfigure. The other parameters are given in Table 1. Notice that the results for the case with only one CB set in the *y*-direction are also shown. Through comparison to an approximate analytical solution, these results can be used as a verification of the numerical implementation. Specifically, by viewing this model as a layered system, we can approximate and using weighted harmonic and arithmetic averages. For example, for the case of *k*_{CB}/*k*_{m} =10^{−3}, we compute =0.029 and =0.96. These values are in close agreement with the numerical results of 0.032 and 0.98, which provides us with some confidence in the numerical computations.

It is evident in Figure 4a (*k*_{CB}/*k*_{m} =10^{−1}) that, as α increases from 5° to 90°, increases from 0.60 to 0.76, whereas decreases from 1.0 to 0.76. Thus, for this value of *k*_{CB}, the impact of the CBs is relatively slight. More impact is apparent for lower values of *k*_{CB}. Specifically, for *k*_{CB}/*k*_{m}=10^{−2}, increases from 0.14 to 0.23, while decreases from 1.0 to 0.23 as α is varied. For *k*_{CB}/*k*_{m} =10^{−3}, increases from 0.016 to 0.028, whereas decreases from 0.90 to 0.028. Note that, for α=90°, , as would be expected. In addition, the results in this limit are close to those for the case with only one CB set in the *y*-direction, which is again consistent with intuition.

The magnitudes and directions of the principal components of the upscaled permeability tensor are also of interest in some applications (recall that the values shown in Fig. 4 are permeability components in the co-ordinate directions and do not fully define the upscaled permeability tensor). Figure 5 shows normalized upscaled permeability tensors using ellipses, with the long- and short-axes representing the major and minor principal permeabilities. These results are computed using area averaging (Equations 5–7). We set *k*_{CB}/*k*_{m} =10^{−2} and consider various intersection angles, α. The magnitudes are normalized by the matrix permeability. As would be expected, for α=5°, the major principal component is nearly aligned with the *y*-direction. As α increases from 5° to 90°, the magnitude of the major principal component decreases from 1.0 to 0.23, and its direction rotates toward the 45° direction (see Fig. 5). The magnitude of the minor principal component increases slightly, from 0.13 to 0.23. For the cases with *k*_{CB}/*k*_{m}=10^{−1} and *k*_{CB}/*k*_{m}=10^{−3}, the trends of the upscaled permeability are similar to those in Figure 5, although the magnitudes of the principle components are significantly different.

### One vertical and one inclined compaction band set

The system considered above was characterized by two sets of vertical CBs, so 2D simulations were appropriate. Previous outcrop observations have shown that vertical (high-angle) CBs and inclined (low-angle or even subhorizontal) CBs can, however, both be present (Aydin *et al.* 2006; Aydin & Ahmadov 2009; Eichhubl *et al.* 2010). Thus, we now consider the impact on flow of a set of inclined CBs intersecting a set of vertical CBs.

The dimensions of the region under study is 20×20×5 m. The vertical CB set for these models was obtained by vertically extending the CB set used for the 2D cases above. This vertical CB set strikes in the *y*-direction. The other (inclined) set was obtained by rotating an identical CB set, with strike in the *x*-direction, on an axis oriented in the *x*-direction that extends through the centre of the model (this axis is defined by the end points *x*=0, *y*=10, *z*=−2.5 and *x*=20, *y*=10, *z*=−2.5). The CB spacing is not changed during these rotations. Figure 6 shows three of the resulting 3D model configurations. Here the thick lines denote CBs, and the rotation (dip angle θ) of the inclined CB set is 0° (upper), 45° (middle) and 90° (lower), respectively. The dip angle, θ, of the inclined CB set is varied from 0° (subparallel to horizontal bedding) to 90° (perpendicular to bedding) in our computations.

The models are discretized using from about 200 000 to about 360 000 cells. In these 3D simulations, the rock matrix is represented by tetrahedral cells and the CBs are represented by (planar) triangular cells. For some cases (particularly at low values of θ), very slight shifts in CB locations are introduced to facilitate meshing. The parameters used in these simulations are given in Table 1. Note that *k*_{CB}/*k*_{m} =10^{−2}.

Figure 7 shows the upscaled permeability components in the *x*-, *y*- and *z*-directions, , and , normalized by the matrix permeability, *k*_{m}, for various dip angles θ of the inclined CB set. We note first that, for θ=90°, the 3D result for and should correspond to the 2D result for α=90° (this 2D configuration is shown in Fig. 3c). In this limit, the 3D computation gives , which agrees with the 2D result for α=90° (Fig. 4b). This provides a degree of verification for the 3D computations. From Figure 7, we see that, as θ increases from 0° to 90°, decreases significantly, increases significantly and stays nearly constant. This is because the strikes of the inclined CBs are in the *x*-direction, which means that the variation of the dip angle of this CB set has a dramatic effect on flow in the *y*- and *z*-directions but only a slight impact on flow in the *x*-direction. These results demonstrate that the flow effects of a vertical CB set plus an inclined CB set can be considerably different from those of two vertical CB sets.

### Joints intersecting compaction bands

As shown in Figure 1c, joints often intersect CBs in the Aztec Sandstone. We now investigate the flow behaviour of such systems. Since we do not have detailed field statistical data on such a system, we combine the field data for the CB set used in the previous examples with joint sets that are constructed using a geologically based statistical method.

We now describe the statistical procedure for generating the joint array. To apply this method, a joint spacing distribution and a joint length distribution are required as inputs. The spacing distribution may be described as evolutionary, as it is dependent on the stage of joint development. Using field, experimental and numerical data, Rives *et al.* (1992) found that the joint spacing distribution is a negative exponential for poorly developed patterns at the initial stage of joint development, log-normal at an intermediate stage and tends towards a normal distribution at a final stage. In this study, we use an interaction process proposed by Rives *et al.* (1992) to generate stochastic joint sets. With this approach, it is assumed that there is a small zone around each joint within which no new joint can form. The small zone (also referred to as the shadow zone in structural analyses) is mainly used to account for the stress relaxation in the vicinity of a joint due to fracturing.

In this work we specify that the shadow zone is within 0.2 m of each joint surface. This value is obtained by multiplying the ratio of fracture spacing to bed thickness (which is assumed to be 0.5 based on Wu & Pollard 1995), with the bed thickness which, in the Aztec Sandstone, is estimated to be 0.4 m (as shown in Fig. 1). In addition, we use a total joint length to control the average joint spacing based on an area method proposed by Wu & Pollard (1995). According to the area method, the average spacing of each joint set is approximately equal to the study area divided by the total length of fractures within that area. The area method is more appropriate than the conventional scanline method for poorly developed joint patterns, and it works equally well for well-developed patterns.

To generate a stochastic joint set, the distribution of joint lengths is also needed. We have a database of the lengths of 120 joints from the Valley of Fire State Park (unpublished data from de Joussineau & Aydin and Davatzas & Aydin). These joint lengths are found to approximately obey a negative exponential distribution, as shown in Figure 8.

The process to generate stochastic joint sets can now be described as follows: (1) generate a random point as the centre of a new joint; (2) check to ensure that the random point is not located in the stress relaxation zone of any other existing joints; (3) sample from the joint length distribution to obtain a (random) joint length; (4) check to ensure that the new joint does not extend into the stress relaxation zone of any existing joint; and (5) repeat the above process until the total joint length reaches an assumed value that is used to control the joint spacing.

Using the above procedure, we generated three joint sets by fixing the total joint lengths to be 200, 600 and 1000 m in the model domain. These correspond to poorly developed, medium-developed and well-developed joint sets, respectively. We combine these joint sets with the CB set used previously, according to the order of their formation; that is, the joints are always younger than the CBs in the Valley of Fire. Figure 9 shows three model configurations, in which the CB set (thick lines) is fixed and oriented in the *y*-direction, and the well-developed joint set (thin lines) is subparallel (upper, α=0°), acute (middle, α=45°) or perpendicular (lower, α=90°) to the CB set. In our computations, we vary α from 0° to 90° and consider all three joint intensities.

Figure 10 shows the upscaled permeability components, normalized by the matrix permeability, for the three joint sets intersecting the CB set at various values of α. For these computations, we take *k*_{CB}/*k*_{m} =10^{−2}. The results for the CB set only are also shown (dashed lines), and it is clear that the presence of the joints impacts the results significantly. We see that, as α increases from 0° to 90°, increases substantially for all three joint sets, whereas decreases. This occurs because the joints act as conduits for fluid to cross the CBs. It is also evident that these effects are magnified for denser joint distributions.

We now repeat these computations using *k*_{CB}/*k*_{m}=10^{−1} and 10^{−3}. Results are shown in Figures 11 and 12. By comparing the results in Figures 10, 11, 12, it is apparent that the magnitude of *k*_{CB}/*k*_{m} has a clear effect on , especially for low α, although the effect on is relatively small. When α is relatively large (e.g. the joints are nearly normal to the CBs), the magnitude of *k*_{CB}/*k*_{m} affects more for the case of poorly developed joints than for medium- or well-developed joints. This is because medium- and well-developed joints provide more conduits for flow to cross CBs, so the magnitude of the CB permeability is less important.

Figure 13 presents the normalized permeability tensors (computed using area averaging) for various combinations of joint and CB sets, for the case *k*_{CB}/*k*_{m} =10^{−2}. In the following, we use and to denote the major and minor principal permeability components, respectively. When the intersection angle, α, is 0° (i.e. joint and CB sets are subparallel), and are in the *y*- and *x*-directions, respectively. For these cases, is about 2–4 times larger than *k*_{m}, whereas is about a quarter of *k*_{m}. Both and increase as the joint density increases. When α is 90° (i.e. joint and CB sets are perpendicular), and are in the *x*- and *y*-directions, respectively. In this case, the principal component in the *x*-direction increases with joint density, while the component in the *y*-direction is always essentially equal to *k*_{m}. As α increases from 0° to 90°, the direction of rotates from the *y*-direction to the *x*-direction. The general trends are similar for other values of *k*_{CB}/*k*_{m}, although the magnitudes of the principal components are different.

To examine the impact of such joint–CB systems on oil recovery, we now present waterflood simulations for the well-developed joint set intersecting the CB set perpendicularly. The rock matrix is assumed to be initially saturated with oil. In the waterflood simulations, we specify an oil viscosity of 2 cP and a water viscosity of 1 cP. Other parameters used in the simulations are given in Tables 1 and 2. We consider waterflood in the *x*-direction, in which water is injected in at the left edge of the model, and oil and water are produced at the right edge of the model, as well as waterflood in the *y*-direction (where flow is from bottom to top). In both cases, we specify fixed pressures at the injection and production edges of the model.

Figure 14 illustrates the water-saturation profiles at 0.25, 0.5 and 0.75 PVI (pore volume injected) for these two simulations. It is apparent that, for the *x*-direction waterflood, joints act as pathways for fluid to cross the CBs. However, because the joints are not interconnected, the injected water contacts the oil in the matrix and provides a reasonably efficient sweep. For waterflood in the *y*-direction, the sweep is very efficient (and is essentially not affected by the joints) because all of the joints are oriented orthogonal to the flow direction.

Figure 15 shows the water cut (defined as the fraction of water in the produced fluid) and the total liquid production rate divided by the imposed Δ*p* (here Δ*p* is the pressure at the left edge minus the pressure at the right edge) for the *x*-direction waterflood. The presence of joints clearly impacts both water cut and total liquid rate, although the level of development of the joints has only a minimal effect on water cut. This is not the case for total liquid rate, which clearly varies for the different joint sets. Of the cases with joints, the system with the well-developed joint set intersecting the CB set is probably the most favourable for waterflood as it gives the highest flow rate and about the same water cut as the other joint–CB cases.

### Small fault overprinting compaction bands

Our final example involves a small fault with one CB set and one joint set, as shown in Figure 1d. Our model of this system is 3D, with the small fault extended vertically by 5 m. Figure 16 shows the structural components considered, where the blue elements represent slip surfaces, between which there is a very thin body of fault rock. As the fault is only about 11 m long and has a slip of a few tens of centimetres, the fault rock is thin (about 3 cm) and discontinuous. Owing to shearing, a set of splay joints (red surfaces) was generated on the right-hand side of the fault forming an asymmetric damage zone. These splay joints also cut through an older set of CBs (green surfaces). We note that it is well established that the faults in the area were formed by the shearing of earlier joints that are subparallel to the pervasive CBs (Flodin & Aydin 2004; Myers & Aydin 2004).

Using the parameters in Table 1, the normalized upscaled permeability components , and are calculated as 0.36, 1.27 and 2.35, respectively. In the *x*-direction, the fault rock and CBs represent barriers for fluid flow, although the splay joints do provide some high-permeability pathways. The slip surfaces and joints are conduits for flow in the *y*- and *z*-directions (especially in the *z*-direction, where portions of the slip surfaces extend over the entire model thickness). Thus, we expect the smallest upscaled permeability to be and the largest upscaled permeability to be , as is observed. Note that these quantities do depend on the size of the upscaled domain.

Both the *x*- and *y*-direction waterflood simulations are also performed for this case. Figure 17 shows the water-saturation profiles for 0.25, 0.5 and 0.75 PVI. In the *x*-direction waterflood, the injected water enters the right part of the domain through discontinuities in the fault rock. Injected water then flows through the splay joints to pass across the CBs. For the *y*-direction waterflood, the slip surfaces and splay joints directly impact the displacement from the very beginning of the simulation (Fig. 17b). The water cut and total liquid production rate are shown in Figure 18. It is apparent that, consistent with Figure 17, water breakthrough occurs earlier for *y*-direction waterflood than for *x*-direction waterflood. The total liquid production rate is higher for the *y*-direction waterflood, consistent with the results for upscaled permeability reported earlier. Note that the imposed Δ*p* is the same for the two cases, although the system lengths differ in the *x*- and *y*-directions.

## Discussion

We have demonstrated that the impact of multiple sets of compaction bands (CBs) on fluid flow and oil recovery can be very different from that of only one set of CBs. The intersection angle between two sets of CBs has a significant effect on the components of the upscaled permeability. Therefore, for many flow simulation applications it is important to incorporate all CB sets into the geological model. In addition, the 3D structure (e.g. the dip angle of inclined CBs) of the CB sets can also impact upscaled permeability. In these situations, the flow problem should not be simplified from 3D to 2D. Previous investigators have reported that the CB permeability may vary from a factor of 10^{−1} to 10^{−4} of that of the surrounding rock matrix (Sternlof 2006; Aydin & Ahmadov 2009; Torabi & Fossen 2009; Sun *et al.* 2011). We found that the CB permeability value can have a significant impact on the upscaled permeability, so it is important that this parameter be accurately quantified.

In our models of joints intersecting CBs, the statistically generated joints were not connected. We showed that, although these joint sets do not themselves interconnect, they still enhance the equivalent permeability normal to the CB set by allowing fluid to cross CBs through high-permeability pathways. If a network of connected joints were to overprint a set of CBs, the impact on upscaled permeability and waterflood performance would be even more substantial.

In this work we considered only one configuration involving a fault. The specific configuration involved the overprinting of a small fault on a set of CBs. Because the fault rock and slip surfaces of the small fault were discontinuous, and because the associated damage zone was small and asymmetric, this fault had less of an impact on flow than would larger faults. The impact of large faults overprinting CBs will be investigated in future work. For such cases, a considerably higher number of associated structural elements will present computational challenges. We note further that, in our simulations, we modelled the slip surfaces of the small fault as fluid conduits with high permeability. However, as shown in Ahmadov *et al.* (2007), if the slip surfaces become slip bands as a result of infilling, they can act to significantly reduce the upscaled permeability.

Our computations were over regions that are roughly the size of a single grid block in a large-scale flow simulation model. Thus, the quantities computed can be used directly in such models. Our results may also be useful for providing appropriate ranges for permeability adjustments during the history matching or data assimilation stage of a field study.

As noted earlier, in all of our finite-volume computations we used two-point flux approximations. This approach can lead to inaccuracy if the grid is significantly nonorthogonal or if the fine-scale permeabilities are anisotropic. To achieve accurate numerical solutions in such cases, the use of multipoint flux approximations, as described by, for example, Aavatsmark *et al.* (1998), will be required.

We note finally that, because CBs and shear bands are similar from the viewpoint of their effect on subsurface flow, our findings for CBs should also apply to shear bands. This would provide our results with broad applicability as shear bands are known to occur widely throughout the world.

## Conclusions

We performed both 2D and 3D flow computations to investigate the impact of the overprinting and intersecting of various types of geological structures on fluid flow in the Aztec Sandstone. The following specific conclusions can be drawn from this study.

The intersection angle of multiple sets of vertical compaction bands (CBs) and the dip angle of inclined CBs have a strong effect on the upscaled permeability tensor. Thus, these features can significantly affect fluid-flow behaviour, such as flow rate (for a given imposed pressure drop) and waterflood performance. The accurate representation of CB patterns and orientation in fluid-flow models is therefore very important.

The intersection angle between joint sets and CB sets strongly influences the principal directions and values of the upscaled permeability tensor. If a set of joints is parallel to a set of CBs, the effect of the joints on the upscaled permeability in the direction normal to the CBs is not significant. By contrast, if a joint set intersects a set of CBs, particularly if it does so at a high angle, the joints provide high-permeability pathways across the CBs that can counteract the impact of the CB set on flow. For the (unconnected) joint sets considered, water cut as a function of PVI is insensitive to the level of development of the joints, although flow rate increases as joint density increases.

For the case of a small fault overprinting a set of CBs, the combined effect of fault rock, slip surfaces, splay joints and CBs acts to reduce equivalent permeability in the fault-normal direction and to increase permeability in the vertical direction. The impact of the fault will, in general, be much greater for larger faults with more and larger size structural elements.

Taken in total, the results in this paper clearly demonstrate the importance of including common types of geological structures with various geometries and flow properties in subsurface flow simulations. Future work should target the consideration of CBs in different orientations forming specific structural domains, and faults of larger size overprinting CBs.

## Acknowledgments

This research has been supported partially by DOE Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences, Geosciences Research Program, grant no. DE - FG03 - 94ER14462, and the Stanford Rock Fracture Project. John Cosgrove and an anonymous reviewer are thanked for their helpful comments. We also thank Jonathan Redfern for inviting and encouraging us to contribute to this volume.

- © The Geological Society of London 2014