Upscaling two-phase flow in naturally fractured reservoirs

Upscaling two-phase flow in naturally fractured reservoirs Stephan K. Matthäi and Hamidreza M. Nick

ABSTRACT Simulation grid blocks of naturally fractured reservoirs contain thousands of fractures with variable flow properties, dimensions, and orientations. This complexity precludes direct incorporation into field-scale models. Macroscopic laws capturing their integral effects on multiphase flow are required. Numerical discrete fracture and matrix simulations show that ensemble relative permeability as a function of water saturation (kri[Sw]), water breakthrough, and cut depend on the fraction of the cross-sectional flux that occurs through the fractures. This fracture-matrix flux ratio (qf/qm) can be quantified by steady-state computation. Here we present a new semianalytical model that uses qf/qm and the fracture-related porosity (ff) to predict kri(Sw) capturing that, shortly after the first oil is recovered, the oil relative permeability (kro) becomes less that that of water (krw), and krw/kro approaches qf/qm as soon as the most conductive fractures become water saturated. To include a capillary-driven fracture-matrix transfer into our model, we introduce the nonconventional parameter Af,w(Sw), the fraction of the fracture-matrix interface area in contact with the injected water for any grid-block average saturation. The Af,w(Sw) is used to scale the capillary transfer modeled with conventional transfer functions and expressed in terms of a rate- and capillary-pressuredependent kro. All predicted parameters can be entered into conventional reservoir simulators. We explain how this is accomplished in both, single- and dual-continua formulations. The predicted grid-block-scale fractional flow (fi[Sw]) is convex with a near-infinite slope at the initial saturation. The upscaled flow equation therefore does not contain an Sw shock but a long leading edge, capturing the progressively widening saturation fronts observed in numerical experiments published previously.

Copyright ©2009. The American Association of Petroleum Geologists. All rights reserved. Manuscript received May 3, 2009; provisional acceptance June 19, 2009; revised manuscript received July 26, 2009; final acceptance August 3, 2009. DOI:10.1306/08030909085

AAPG Bulletin, v. 93, no. 11 (November 2009), pp. 1621–1632

1621

AUTHORS Stephan K. Matthäi Reservoir Engineering, Department of Mineral Resources and Petroleum Engineering, Max-Tendler Strasse 4, Montan University of Leoben, Leoben 8700, Austria; [email protected] Stephan K. Matthäi is the chair of reservoir engineering at the University of Leoben, Austria. Before that, he was a senior lecturer of computational hydrodynamics at Imperial College London, and a postdoctoral researcher at Eidgenössische Technische Hochschule Zurich, Switzerland; Stanford University; and Cornell University. His Ph.D. is from the Research School of Earth Sciences, Australian National University. He holds a diploma degree from Eberhard Karl’s University, Tübingen, Germany. Hamidreza M. Nick Center of Petroleum Studies, Department of Earth Science and Engineering, Imperial College London, London, United Kingdom; [email protected] Hamidreza Maghami-Nick is currently a Ph.D. student in the Center of Petroleum Studies at Imperial College London. He holds M.Sc. degrees from the Utrecht University, Netherlands, and from the K.N.T. University of Technology, Iran. His research interests range from the development of numerical methods to upscaling solute transport in fractured porous media.

ACKNOWLEDGEMENTS This work was generously supported by the sponsors of the United Kingdom Industry Technology Facilitator and the Technology Strategy Board consortium “Improved Simulation of Fractured and Faulted Reservoirs,” Phase 2. We thank O. Gosselin and an anonymous reviewer for their constructive comments, which helped improve the original manuscript. The AAPG Editor thanks the following reviewers for their work on this paper: Olivier R. Gosselin and an anonymous reviewer.

INTRODUCTION The inability of conventional simulation tools to forecast production from naturally fractured reservoirs (NFRs) has sparked a significant volume of research into alternative simulation methods. One of these is discrete fracture simulation on unstructured grids (e.g., Kiraly and Morel, 1976; Sudicky and McLaren, 1992). New discrete fracture network (fracture-only) models (DFN) (e.g., Dershowitz and Miller, 1995; Unsal et al., in press) and discrete fracture and matrix (DFM) models were developed (e.g., Kim and Deo, 2000; Juanes et al., 2002; Geiger et al., 2004; Matthäi et al., 2005) and are used increasingly to compute the ensemble permeability of fractured rocks (e.g., Matthäi and Belayneh, 2004; Matthäi et al., 2007b). The main unknowns to be addressed by field-scale simulations are as follows. • Productivity of wells and sustainability of flow rate • In the case of injection, time to injectant breakthrough (see Figure 1a where production is halved following water breakthrough) • Recovery (shape of the time-integrated relative permeability curves, Figure 1b) • Recovery mechanism(s) controlling rate dependency • Residual saturations and their distribution However, all fractures cannot be included into fieldscale models because there may be billions of them in each cubic kilometer of reservoir rock. Upscaling is therefore a prerequisite for field-scale simulation. At the same time, kilometer-scale fractures like in corridors require an explicit deterministic representation; otherwise, the field-scale model will not be predictive (e.g., Singh et al., 2008). Existing upscaling approaches focus on permeability (k) (e.g., Oda, 1985; Bogdanov et al., 2007; KarimiFard et al., 2006) and/or are directed at deriving paramFigure 1. Uncertainties associated with the production of naturally fractured reservoirs. (a) Typical water breakthrough curve and (b) uncertainty in the reservoirscale relative permeability of oil as illustrated by the gray-shaded area. Terms are defined in Table 1.

1622

Upscaling Two-Phase Flow

eters for dual-continua simulations (e.g., Sarda et al., 2002). To conduct multiphase flow simulations on NFR field-scale models, scale transformations of the capillary pressure (pc), water saturation (Sw), relationship, and relative permeability (kri[Sw]) are required because these originate from the interpretation of multiphase flow experiments on the core-plug scale where capillarity continuity can be assumed and where conventional kri(Sw) models have a physical meaning. However, because such small-scale variables also serve as input to conventional simulators, we advise to stay within this well-established framework so that upscaled properties can be input into these tools. This leaves two options for the transformation: (1) single-porosity upscaling and (2) dual-continua upscaling. In this article, we will show how upscaling can be accomplished with the help of the fracture-matrix flux ratio (qf/qm), and the fraction of the total fracture-matrix interface area in contact with water as a function of gridblock average saturation (Af,w[Sw]). We will show that the latter parameter can also be expressed by a gridblock saturation-dependent shape factor, S = (AfAf,w [Sw]/V)(1/L), where V is the specific volume and L the mean matrix block diameter (see Table 1). Our two-phase flow upscaling method is informed by numerical DFM simulations of viscous two-phase flow as will occur during water flooding of a fractured reservoir (Matthäi et al., 2007b). This analysis is restricted to flow regimes where gravitational forces are insignificant and capillary fracture-matrix transfer is not modeled. However, capillary transfer is still considered by our new upscaling formulation: this is accomplished using standard transfer functions applied to modify the relative permeability of the nonwetting phase as a function of capillary forces and flow velocity. Capillary-driven fracture-matrix transfer is therefore accounted for through modification of the grid-block-scale sweep efficiency. This leads to the frequently observed rate

Table 1. Geometrical and Flow Properties Derived by Analysis of DFM Models Notation

Units

fm Vm ff Vf ffX Af

m3 m3 X m2

D k k vt v qf/qm lt kri(Sw)

m m2 m2 m s−1 m s−1

Swc Snr Swe T(…) Af,w(Sw)

X X X X X

Swf(Sw) Subscript Subscript Subscript Subscript Subscript BC CCI

X eff f m o w

Name

Discretization

Pore volume of rock matrix Void space volume in fractures Fraction of total f caused by fractures Specific fracture-matrix interface area (equals to P32 in DFN modeling literature) Mean or mode of block diameter Permeability tensor Permeability in the direction of flow Darcy velocity of fluid mixture Scalar velocity magnitude Fracture-matrix flux ratio Total mobility Relative permeability of phase i as a function of average water saturation Connate or irreducible water saturation Irreducible saturation of nonwetting phase Effective water saturation Fracture-matrix transfer function Fraction of fracture-matrix interface area in contact with water as a function of average water saturation Fracture saturation as a function of average water saturation Effective With regard to fracture With regard to rock matrix With regard to oil With regard to water Superscript denoting parameters calculated with the Brooks-Corey relative permeability model Superscript referring to countercurrent imbibition

dependency of large-scale relative permeability, manifesting itself in the recovery curves; the time to water breakthrough; and the water cut. Importantly, if upscaling could be accomplished solely by the adjustment of parameter values, large-scale system behavior would be readily predictable. However, game theory and complex systems science have established that complex systems like NFRs exhibit emergent behavior. This cannot be predicted from component behaviors studied process by process (cf. Holland, 1997). The challenge in upscaling multiphase flow in NFRs therefore lies in capturing this emergent behavior. This article is structured as follows. First, all steps of the DFM modeling process are described establishing an analysis strategy for such experiments. Key observations on the model behavior are reviewed and are

Piecewise Piecewise Piecewise Piecewise

constant constant constant constant

Piecewise constant Piecewise constant Piecewise constant Piecewise constant Piecewise constant NA Piecewise constant Piecewise constant Piecewise constant Piecewise constant Piecewise constant Piecewise constant NA Piecewise constant

related to experiments described in the literature. Second, our new model for the grid-block-scale relative permeability is presented. This model is then enriched with transfer functions considering capillary fracturematrix transfer. Third, we show how the derived parameters can be used for single- and dual-porosity upscaling. Finally, we apply our new model on the larger scale and observe the ensuing displacement patterns and volumeaveraged behavior.

METHODS In this section, we set out the workflow and analytical methods we have developed to measure fracture-matrix Matthäi and Nick

1623

ensemble single and multiphase flow properties of DFM models built from outcrop analogs or fracture statistical data obtained from well logs and core samples. This workflow has the following steps. 1. Fracture modeling: collect data and build a boundary representation of a grid-block-scale model of fractured rock using a computer-aided design (CAD) system and smooth parametric curves and surfaces (NURBS) to capture the geometry of fracture surfaces and bedding planes (e.g., Matthäi and Belayneh, 2004; Belayneh et al., 2006). 2. Unstructured spatially adaptive finite-element meshing of the CAD model (see Paluszny et al., 2007). 3. Geomechanical modeling of fracture aperture distributions (described in the companion article and by Paluszny and Matthäi, 2009). 4. Steady-state analysis of fracture porosity, specific fracture-matrix interface area, effective permeability, and fracture-matrix flux ratio and their variability with orientation of the far-field stress (Matthäi et al., 2007b). 5. Forward simulation of imbibition and drainage to determine the total mobility and dynamic relative permeability-average saturation relationships (Matthäi et al., 2005, 2007a). The models are simple. We use standard relative permeability-saturation relationships for the rock matrix such as the Brooks-Corey model (Brooks and Corey, 1964) or experimentally determined ones (cf., Belayneh et al., 2006). For the fractures, we apply linear relative permeability-saturation relationships except for those with apertures less than 0.5 mm (0.01 in.) (approaching the diameter of the largest pores in the host rock); here we use the BrooksCorey model with a very small exponent, as discussed by Helmig (1997). 6. Curve fitting and semianalytical representation of the determined ensemble properties with the new upscaled relative permeability model. 7. Postprocessing of parameters for use in either singleor dual-porosity simulation. 8. Inter- or extrapolation of upscaled properties across reservoir layers. This can be conducted with the standard reservoir characterization and modeling tools. For this purpose, correlations between fracture properties and geological attributes of the layers should be analyzed. A systematic study of how the orientation of the principal axes of local qf/qm tensors depends on the dip/dip direction of the fractured layer and the in-situ stress would be most helpful here. 1624


The main scope of this article is the upscaling method (steps 6 and 7), but the analysis of the underpinning numerical DFM experiments (steps 4 and 5) is also explained in some detail. The remainder of the workflow is explained elsewhere and is referenced when possible.

Property Extraction from Numerical DFM Experiments Paluszny et al. (2007) and Matthäi et al. (2007b) laid out how hectare-meter-scale DFM models can be built and analyzed for their effective permeability (keff), fracture-matrix flux ratio (qf /qm), total mobility (lt), and ensemble relative permeability (kri[Sw]). The numerical determination of kri is achieved with a hybrid finiteelement node-centered finite-volume simulation technique (Paluszny et al., 2007), facilitating two-phase flow simulations with complex fracture geometries based on outcrop analogs (Belayneh et al., 2006) or geostatistical fracture models (Belayneh et al., 2009, this issue). Whether a stochastic model is sufficient depends on the degree of self-organization caused by fracture interactions during growth (cf. Wu and Pollard, 1995). If these were significant, geomechanical models are required to predict fracture connectivity correctly (e.g., Renshaw, 1997, Paluszny, 2008). We recommend that geomechanical modeling be used at least to compute the fracture aperture and permeability distributions considering the in-situ stress state in the reservoir. Mounting evidence is observed from field studies and numerical experiments that the in-situ stress exerts a decisive influence on the transmissivity of fractures (Barton et al., 1995; Zhang and Sanderson, 1996). The properties extracted by the analysis formalized in our modeling workflow are summarized in Table 1. Effective Permeability and Fracture-Matrix Flux Ratio (Workflow Step 4) To compute keff and qf /qm, the single-phase steadystate fluid pressure equation

keff 5p ¼0 5 m

ð1Þ

is solved for uniform fluid viscosity (m) and fluid pressure (p) values assigned to opposite model boundaries. We accomplish this with a linear finite-element method applied

to the fully discretized fracture and matrix model using node-centered finite volumes to integrate the fluid flux across model boundaries. Inclusion of the rock matrix is very important because in reservoir rocks, almost by definition, significant viscous flow through the rock matrix occurs (Matthäi and Belayneh, 2004). The computation of qf/qm involves the calculation of the matrix volume flux in the absence of fractures, qm = A(hkmi/m)∇p, ensuing from a specific far-field fluid-pressure gradient, ∇p = (pb1 − pb2)/L, assigned to the DFM model with length L and cross-sectional area A (see Matthäi and Belayneh, 2004). Variable hkmi is the average matrix permeability. The box-shaped model should be aligned with the principal axes of the permeability tensor. The qf/qm ratio can then be determined from qm ¼

fm Vm Akm rp ff Vf þ fm Vm m

ð2Þ

and the total flux (qt) qf ¼ ðqt qm Þ=qm qm

ð3Þ

Ensemble Relative Permeability Conventional steady-state relative permeability measurements are conducted under fully satiated conditions; i.e., both fluid phases form a percolating cluster across the sample at all times (Loomis and Crowell, 1962). Here we use a different unsteady approach to determine kri(Sw), starting in the single-phase regime, except for when irreducible matrix saturations exist. We follow this route because we are also interested in the time to displacing-phase breakthrough (Matthäi et al., 2007b). This methodology is also known as dynamic relative permeability measurement (e.g., Dullien, 2004; Schembre and Kovscek, 2006). Relative permeability values are calculated from the incremental saturation change for each fluid volume injected (Qt) as measured by the finite-volume surface integration of the piecewise constant finite-element total velocity (vt) over the inflow boundary held at constant saturation. Both ensue in response to the applied far-field fluid pressure gradient. This leads to the average boundary-normal (scalar) velocities nt ¼

Qt ADt

nt ¼ no þ where fm Vm, and ff Vf are the matrix and fracture pore volumes, respectively. Both keff and qf /qm are tensor properties with offdiagonal terms when their principal axes are not aligned with the coordinate system.

Ensemble Total Mobility In simulations with constant pressure boundary conditions and a known model keff as computed with the procedure described above, the total mobility (lt) as a function of the average water saturation (Sw) is inferred from qt ¼ lt Akeff rp; lt ¼

qt Akeff lp

ð4Þ

using incremental measurements of the volumetric flow rate, qt (m3 s−1) per unit cross-sectional area. In the singlephase end-member cases, lt is simply 1/mi. We stress that the mesoscopic fluid pressure (p) in equation 4 is that of the continuous wetting phase. Capillary-driven flow is not modeled by these numerical experiments, as will be addressed later. In their analysis of DFM models, Matthäi et al. (2007b) identified a nonlinear, consistently low lt as a key characteristic of a fractured rock.

DSw ðfm Vm þ ff Vf Þ ADt

ð5Þ

From the extended Darcy’s law and the fluid viscosities (mi), we find kri ∣vw ∣mw ; keff rp ∣vo ∣mo kro ðSw Þ ¼ keff rp

krw ðSw Þ ¼

ð6Þ

This is only one of the different possibilities concerning how relative permeability-saturation relationships can be extracted from the numerical experiments. The key challenge associated with this procedure is that fracture porosity tends to be very low (0.001–1%). However, if the fractures dominate the flow (qf/qm > 1), then their saturation state controls kri(Sw). Thus, kri(Sw) is extremely sensitive to small saturation changes. Therefore, experimental kri(Sw) data may appear erratic (Figure 2) and are difficult to fit, as has already been addressed by Matthäi et al. (2007b). At this stage of the analysis, a decision needs to be made on whether the results are upscaled into a single- or a dual-continuum model. If a single-continuum model were used, the analysis would proceed with the determination of the residual saturations in imbibition and drainage runs. Matthäi and Nick

1625

Figure 2. Experimentally derived ensemble relative permeability curves. (a) Single-flow aligned fracture in porous matrix (redrawn from Haghighi et al., 1994; reprinted with permission from Elsevier); (b) our own DFM results for a fracture geometry built with data from Clair field, Shetland islands, United Kingdom; and (c) multiple fractures in a porous matrix (from Maloney and Doggett, 1997; public domain). Note the rapid crossover as water forms a continuous path across the sample. Terms are defined in Table 1.

For upscaling into dual-continua models, we recommend monitoring water saturation in the fractures and the rock matrix, separately, obtaining Swf, and Swm. Because kri(Sw) is already known for the rock matrix, an ensemble kri(Sw) can be found as a qf/qm-weighted average of model kri(Swf) and matrix kri(Sw). The function kri(Swf) is much less erratic than its single-porosity counterpart but is constant over a wide range of model saturations. This and other relative permeability issues are further discussed in a subsection in conjunction with the introduction of our new upscaled relative permeability model. The Af,w(Sw) is the fraction of the fracture-matrix interface area in contact with the injected fluid (Figure 3) as a function of Sw and normalized by the total fracturematrix interface area (Af). The Af,w(Sw) is measured simultaneously for all fractures, using Sw values interpolated from the finite-element nodes (node-centered finite volumes) to their barycenters. Where these exceed a certain threshold value, the interfacial area of the finite-element faces situated at the fracture-matrix interface is integrated. If lower dimensional (surface) elements are used to represent the fractures (cf. Juanes et al., 2002), their area is used directly but is counted twice because each fracture segment has matrix interfaces on either side of 1626


the fracture void. Thresholding of fracture Sw is applied because the saturation front in the fracture can be highly dispersive; indeed, when a linear fracture-relative permeability model is used, the measurement of Af,w(Sw) requires the definition of a specific threshold value distinct from the initial equilibrium fracture saturation. Considering the symmetry of the front-normal saturation profile, we choose a value of half of the inlet saturation. Thresholding and saturation interpolation introduce an error into the Af,w(Sw) estimate, which scales with the mesh resolution. Hence, fracture meshes have to be sufficiently refined to obtain meaningful measurements. Fracture Saturations If upscaling into a dual-continuum model is the goal, fracture and matrix saturations can be recorded, separately. In this case, all properties should be recorded as a function of the fracture saturation (Swf) because of the greater sensitivity of the model behavior to this parameter. This is further discussed in a following section. End-Point Saturations The length of the imbibition and drainage runs is controlled either by fixing the maximum number of pore

Figure 3. The FRACS500 quarter of the five-spot problem with a water injector on the left vertical edge and a producer on the right. (a) Blue to green fractures indicate Af,w(Sw) after an overall saturation change of 0.2. Sharp fracture saturation fronts indicate that (b) the Af,w(Sw) estimate as a function of the model average Sw obtained by thresholding is accurate. Terms are defined in Table 1.

volumes injected or by choosing a maximum span of actual time allowed at the flow rate of interest (e.g., 20 yr if this is the duration of a supply contract serviced by the operator). The ensemble residual saturation of the nonwetting phase (Sor) and the connate (irreducible) water saturation (Swc) are measured at the end of such runs. They are bounded by the original values assigned to the rock matrix and the fractures, respectively, and can be used to calculate an ensemble effective water saturation (Swe) Swe ¼

Sw Swc 1 Swc Sor

ð7Þ

so that kri(Sw) can be expressed in terms of kri(Swe). Importantly, visual examination of the saturation distribution model at the end of a run reveals where the displaced phase gets trapped, providing clues as to which enhanced oil recovery method carries the highest potential for recovering this oil.

UPSCALING METHOD The experiments of Matthäi et al. (2007b) and Belayneh et al. (2009, this issue) showed that grid-block relative permeabilities are very sensitive to volume-averaged saturation. Matthäi et al. (2007b) speculated that this behavior could be fitted with a new grid-block-scale relative permeability model. For this purpose, they employed the three parameters qf/qm, Af,w(Sw), and the fraction of the void space caused by the presence of the fractures (ff). However, capillary-driven fracture-matrix transfer is not accounted for by this model. Here we propose to upscale this important process in terms of an extra contribution to the ensemble relative permeability of the nonwetting phase enhancing sweep efficiency. This kro(Sw) contribution has to be scaled by the average Darcy velocity (vt) in the upscaled (homogenized) grid block because it depends on the local balance between capillary and viscous forces. We suggest quantifying the

Matthäi and Nick

1627

capillary contribution to kro(Sw) using two additional terms: the “rate of countercurrent imbibition (T[Sw, Sw,i]),” which is dependent on the fracture-matrix capillary pressure differential, initial grid-block saturation (Swi), and the current average grid-block saturation (Sw). This rate can be calculated with any conventional transfer function of choice. However, the choice of T(Sw,Sw,i) should be consistent with the porous medium model, which is used to calculate the two-phase flow properties of the rock matrix. Thus, if the BrooksCorey kri is used for the rock matrix, it should also be used inside T(Sw,Sw,i). In this context, it is important that the analysis by Zimmermann et al. (1996) shows that, for early times, T(Sw,Sw,i) only depends on Af/V, but for late pseudo-steady-state transfer, a “representative rock-matrix-block diameter” (D) must be specified as well. As in dual-continua models, T(Sw,Sw,i) itself should use (Af/V) (1/D) as a scaling (shape) factor. If fracture spacing is used as a proxy for D, the spacing distribution should be rigorously established, and its mode instead of the mean should be used to quantify D if the distribution is lognormal. Other aspects of transfer functions are reviewed by Abushaikha and Gosselin (2008). The Af,w(Sw) is used to scale T(Sw,Sw,i) calculated for each grid block. The following equations, presented here for the first time, summarize our new upscaling model for viscous and capillary two-phase flow in fractured porous media: krw ðSw Þ ¼

erfðSw =ff ÞSw ð1 Sw Þðqm =qf Þ2 þ Sw

kro ðSw Þ ¼ max 1 krw ½Sw ;

erf½Sw =ff ½1 Sw ½1 Sw þ Sw ½qf =qm 2 !

þ kCCI ro ½Sw kCCI ro ðSw Þ ¼ Af;w ðSw ÞTðSw ; Sw;i Þ=ð1 þ vt qf =qm Þ

ð8Þ

where erf( ) is the standard error function and superscript CCI stands for capillary transfer of oil into fastflowing fractures by countercurrent imbibition. Because early-time capillary fracture-matrix transfer can be very ðSw Þ fast, application of the max( ) function to kupscaled ro guarantees that the volume of nonwetting phase that leaves the grid block is limited by the total flow that occurs through it in any time interval. This condition is not required in the dual-porosity upscaling as further described later, but only if transfer functions are used, which consider fracture saturation. 1628


To achieve a smooth transition between fractured and unfractured rock, we average the upscaled relative (Sw ) from equation 8 with that of permeability kupscaled ri the porous matrix (kBC rmi ½Sw where m refers to the matrix permeability) as calculated here with the Brooks-Corey model.

ðSw Þ ¼ kensemble ri

kupscaled ðSw Þ þ ðqm =qf Þ2 kBC ri rmi ðSw Þ 1 þ ðqm =qf Þ2

ð9Þ

In this averaging, we use the fracture-matrix flux ratio as a weighting factor. This completes the single-porosity upscaling formulation. The new formulation can also be used to upscale two-phase flow using a dual-porosity framework. Now, keff and the upscaled ensemble relative permeability are used for the flowing (fracture) domain. The key difference from the single-porosity formulation presented in equations 8 and 9 is that the dual-continuum model permits keeping track of the water saturation in the fractures (Swf) and the rock matrix (Swm), separately, so they can be used in the transfer functions. In this framework, upscaled (Sw ), whereas Af,w(Sw) is kCCI ro is dropped from kro retained to scale the transfer computed with T(Sw,Sw,i). The shape factor for the transfer function is computed from Af/V (the fracture-matrix interface area per unit volume of fractured rock = P32 of DFN modeling) and the block diameter D, determined by analysis of the fracture data. The dual-porosity transfer term between fractures and rock matrix thus becomes T ¼ Af;w ðSw Þ

Af 1 TðSw;f ; Sw;m Þ=ð1 þ nt qf =qm Þ ð10Þ VD

Note the consideration of the flow-velocity magnitude as in the single-porosity formulation. Figure 4 shows how the fracture-matrix interface area is distributed between the larger and the smaller fractures of the model FRACS2000 of Matthäi et al. (2007b). This model has a power-law diameter-frequency distribution. This model and its transient saturation distribution (Matthäi et al., 2007b, their figure 26) illustrate that fracture saturation will vary greatly among the fractures and that a dual-continuum model will just be able to keep track of an average of this distribution. For this reason, the same scaling of capillary fracture-matrix transfer needs to be applied as already conducted using Af,w(Sw) in the single-porosity upscaling.

Figure 4. Diameter-frequency distribution of the fractures in the model FRACS2000 of Matthäi et al. (2007b). The cumulative (stippled) Af curve shows that most of the interface area is associated with small fractures. However, the velocity histogram (Matthäi et al., 2007b, their figure 23) shows that the flow velocity in these fractures is scarcely any faster than that of the rock matrix. Also, note that the largest fractures with the highest flow velocity contribute little to Af. Terms are defined in Table 1. Although we have some Af,w(Sw) data from DFM simulation, we have not yet identified a functional form of this parameter. For large fracture aperture variations and small mobility ratios, Af,w(Sw) appears to be an exponential of Sw. In the absence of fluid viscosity contrasts and moderate aperture variations, it resembles the cumulative curve in Figure 4 to a greater extent. More sensitivity analyses are required to constrain this parameter.

RESULTS Figure 5 shows kensemble (Sw ) curves obtained for difri ferent fracture qf/qm ratios via single-porosity upscaling (equations 8, 9). These results indicate that as q f /qm decreases, so does the effect of the fractures on (Sw ). kensemble ri

These upscaled relative permeability curves imply a fractional flow function whose derivative no longer switches sign with increasing saturation: it increases monotonically. Observed and predicted grid-block fractional flow is highly dispersive. This is a consequence of the volume averaging intrinsic to any relative permeability formulation but is performed here on a much larger scale than is amenable to laboratory analysis. Fractional flows calculated from the upscaled relative permeabilities (Figure 6) no longer have an inflection point at the water saturation that corresponds to the characteristic shock in the solution of the multiphase flow equation (e.g., Dake, 1978). Instead, the slope of the fractional flow function is maximized near the irreducible saturation and decreases monotonically above this value. This gives rise to the highly dispersive behavior, i.e., the saturation front with a long leading edge. Front-normal cross-sectional averaging of saturation

Figure 5. Upscaled grid-block-scale ensemble relative permeability-saturation relationships for different fracture-matrix flux ratios. Ensemble krw is shown in red, kro in green, kro including CCI in yellow, and sum of the kri in blue. Terms are defined in Table 1. Matthäi and Nick

1629

Figure 6. Upscaled ensemble fractional flow as a function of water saturation as calculated from the ensemble relative permeability-saturation relationships. The stippled line shows the shape of a porous media fractional flow function for comparison. Terms are defined in Table 1.

in DFM models confirms this; this behavior is well matched by predictions using our new model (see Figure 7).

DISCUSSION The rate of countercurrent imbibition is well constrained for the case of countercurrent flow across a material interface (van Duijn and de Neef, 1998) when saturations

Figure 7. Comparison of frontnormal, cut-plane-averaged water saturation and prediction made for the model FRACS2000 of Matthäi et al. (2007b) 162 days after the onset of water injection. Note the long-leading edge that is characteristic of the front. Terms are defined in Table 1.

1630


on either side only reflect this process. In practice, however, viscous, gravitational, and capillary flows will simultaneously affect interface-normal saturation gradients such that this treatment is no longer valid. In a real fractured reservoir, countercurrent imbibition will occur only in those fractures that are actually imbibed by water. Our results indicate that until several pore volumes have been injected, this will only be a small subset of the total number of fractures. The number of fractures that are infiltrated will also depend on the flow rate. The CCI will defend fractures against invasion of the injected fluid so that viscous flow is affected by it, especially in very small fractures, which exhibit nonlinear relative permeability saturation relationships. In the presence of viscosity contrast between the two fluids, CCI might enhance fingering again affecting Af,w(Sw). This dependence is not yet accounted for appropriately by the upscaling formulation and should be subject to future research. The practical implication of this is that for the moment DFM experiments have to be conducted for the specific fluid viscosity contrast and flow rates observed in the reservoir of interest. How are the upscaled parameters distributed spatially on a field-scale reservoir model? Although this is not addressed in this article, the answer to this question should be similar to the handling of special core analysis data. Once a spectrum of parameters has been obtained representative of the different fracture and stress patterns and flow velocities expected in the reservoir of interest, these should be extrapolated onto the layers of the NFR. For the anisotropic fracture-matrix flux ratio, one could use potential relationships between fracture orientation and layering to orient the

tensors. Alternatively, stochastic methods or techniques like sequential indicator simulation could be used to distribute upscaled parameters across models. More research is needed to develop accurate and efficient procedures for this.

CONCLUSIONS This article presents a workflow and semianalytical model for the upscaling of multiphase flow properties of NFRs from the meter to the grid-block scale. The goal is to develop a method to compute multiphase flow properties from field data for use in conventional reservoir simulation where subseismic fractures cannot be represented explicitly. The upscaling method presented is underpinned by numerical simulations of two-phase flow through three-dimensional DFM models built from field data. These show that grid-block-scale relative permeability, time to water breakthrough, and water cut depend primarily on the fracture-to-rock matrix flux ratio (qf/qm), i.e., the fraction of the total crosssectional flux that occurs through the fractures. This tensor property is quantified by steady-state calculations of the partitioning of the Darcy flux between fractures and rock matrix. Using qf/qm and the fracturerelated porosity ff, the ensemble relative permeability is closely approximated by our new ensemble relative permeability model. Capillary-driven fracture-matrix transfer is included into the upscaling by the introduction of the nonconventional parameter Af,w(Sw), the fraction of the fracture-matrix interface that is in contact with water as a function of grid-block average saturation. This parameter is used to scale capillary transfer as modeled with conventional transfer functions. Our semianalytical model represents this transfer as a rate- and capillary-pressure-dependent kro and grid-block-scale sweep efficiency. We have also shown how the upscaled parameters can be used in both single- and dual-continua formulations for reservoir simulation on the field scale.

REFERENCES CITED Abushaikha, A., and O. Gosselin, 2008, Matrix-fracture transfer function in dual-medium flow simulation: Review, comparison, and validation: Europec/European Association of Geoscientists and Engineers Conference and Exhibition, June 9–12, 2008, Rome, Italy, Society of Petroleum Engineers Paper No. 113890. Barton, C., D. Mos, and M. Zoback, 1995, Fluid flow along potentially active faults in crystalline rock: Geology, v. 23, no. 8, p. 683–686,

doi:10.1130/0091-7613(1995)0232.3 .CO;2. Belayneh, M., S. Geiger, and S. K. Matthäi, 2006, Numerical simulation of water injection into layered fractured carbonate reservoir analogs: AAPG Bulletin, v. 90, no. 10, p. 1473–1493, doi:10.1306/05090605153. Belayneh, M., S. K. Matthäi, M. J. Blunt, and S. Rogers, 2009, Comparison between deterministic with stochastic fracture models in water-flooding numerical simulations: AAPG Bulletin, v. 93, p. 1633–1648 Bogdanov, I. I., V. V. Mourzenko, J.-F. Thovert, and P. M. Adler, 2007, Effective permeability of fractured porous media with power-law distribution of fracture sizes: Physical Reviews E, v. 76, no. 3, 15 p., doi:10.1103/PhysRevE.76.036309. Brooks, R. H., and A. T. Corey, 1964, Hydraulic properties of porous media: Colorado State University hydrology paper 3, Fort Collins, Colorado State University Press. Dake, L. P., 2001, 1978, Fundamentals of reservoir engineering, 19th ed.: Developments in petroleum science 8: Amsterdam, Elsevier, 443 p. Dershowitz, B., and I. Miller, 1995, Dual porosity fracture flow and transport: Geophysical Research Letters, v. 22, no. 11, p. 1441– 1444. Dullien, F. A. L., 2004, Porous media: Fluid transport and pore structure, 2d ed.: San Diego, Academic Press, 574 p. Geiger, S., S. Roberts, S. K. Matthäi, C. Zoppou, and A. Burri, 2004, Combining finite element and finite volume methods for efficient multi-phase flow simulation in highly heterogeneous and geometrically complex porous media: Geofluids, v. 4, p. 284– 299, doi:10.1111/j.1468-8123.2004.00093.x. Haghighi, M., B. Xu, and Y. C. Yortsos, 1994, Visualization and simulation of immiscible displacement in fractured systems using micromodels: I. Drainage: Journal of Colloid and Interface Science, v. 166, no. 1, p. 168–179, doi:10.1006/jcis.1994.1283. Helmig, R., 1997, Multiphase flow and transport processes in the subsurface: A contribution to the modeling of hydrosystems: Berlin, Springer-Verlag, 279 p. Holland, J. H., 1997, Emergence: Reading, Massachusetts, Helix Books, Addison Wesley, 258 p. Juanes, R., J. Samper, and J. Molinero, 2002, A general and efficient formulation of fractures and boundary conditions in the finite element method: International Journal for Numerical Methods in Engineering, v. 54, p. 1751–1774, doi:10.1002/nme.491. Karami-Fard, M., B. Gong, and L. J. Durlofsky, 2006, Generation of coarse-scale continuum flow models from detailed fracture characterizations: Water Resources Research, v. 42, W10423, doi:10.1029/2006WR005015. Kim, J. G., and M. D. Deo, 2000, Finite element, discrete-fracture model for multiphase flow in porous media: American Institute of Chemical Engineers Journal, v. 46, no. 6, p. 1120– 1130. Kiraly, L., and G. Morel, 1976, Remarques sur l’hydrogramme des sources karstiques simule par modeles mathematiques: Bulletin du Centre d’Hydrogeologie 1, Institut de Geologie, Universite de Neuchatel, Switzerland, p. 37–60. Loomis, A. G., and D. C. Crowell, 1962, Relative permeability studies: Gas-oil and water-oil systems: Washington D.C., U.S. Government Printing Office, U.S. Bureau of Mines Bulletin 599, 39 p. Maloney, D., and K. Doggett, 1997, Multiphase flow in fractures: BDMP Technologies, Department of Energy Contract DEAC22-94PC91008, Report SCA-9730, 11 p. Matthäi, S. K., and M. Belayneh, 2004, Fluid flow partitioning between fractures and a permeable rock matrix: Geophysical Research Letters, v. 31, no. 7, L07602, doi:10.1029/2003GL019027. Matthäi, S. K., A. A. Mezentsev, and C. C. Pain, 2005, A higher-order TVD transport method for hybrid meshes on complex geological

Matthäi and Nick

1631

geometry: International Journal for Numerical Methods in Fluids, v. 47, p. 1181–1187. Matthäi, S. K., et al., 2007a, Numerical simulation of multi-phase fluid flow in structurally complex reservoirs, in S. J. Jolley, D. Barr, J. J. Walsh, and R. J. Knipe, eds., Structurally complex reservoirs: Geological Society (London) Special Publication 292, p. 405–429. Matthäi, S. K., A. Mezentsev, and M. Belayneh, 2007b, Finiteelement node-centered finite-volume experiments with fractured rock represented by unstructured hybrid element meshes: Society of Petroleum Engineering Journal of Reservoir Evaluation and Engineering, v. 10, no. 6, p. 740–756. Oda, M., 1985, Permeability tensor for discontinuous rock masses: Geotechnique, v. 35, p. 483–495. Paluszny, A., 2008, Numerical modeling of fracture propagation and its implications for fluid flow: Ph.D. thesis, Department of Earth Sciences and Engineering, Imperial College London, London, United Kingdom, 207 p. Paluszny, A., and S. K. Matthäi, 2009, Numerical modeling of discrete multi-crack growth applied to pattern formation in geological brittle media: International Journal of Solids and Structures, v. 46, no. 18–19, p. 3383–3397, doi:10.1016/j.ijsolstr .2009.05.007. Paluszny, A., S. K. Matthäi, and M. Hohmeyer, 2007, Hybrid finiteelement finite-volume discretization of complex geologic structures and a new simulation workflow demonstrated on fractured rocks: Geofluids, v. 7, p. 186–208, doi:10.1111/j.14688123.2007.00180.x. Renshaw, C., 1997, Mechanical controls on the spatial density of opening-mode fracture networks: Geology, v. 25, no. 10, p. 923–926, doi:10.1130/0091-7613(1997)0252.3.CO;2. Sarda, S., L. Jeannin, R. Basquet, and B. Bourbiaux, 2002, Hydraulic

1632


characterization of fractured reservoirs: Simulation on discrete fracture models: Society of Petroleum Engineers Journal of Reservoir Evaluation and Engineering, v. 4, p. 154–162. Schembre, J. M., and A. R. Kovscek, 2006, Estimation of dynamic relative permeability and capillary pressure from counter-current imbibition experiments: Transport in Porous Media, v. 65, p. 31–51. Singh, S. K., H. Abu-Habbiel, B. Khan, M. Akbar, A. Etchecopar, and B. Montaron, 2008, Mapping fracture corridors in naturally fractured reservoirs: An example from Middle East carbonates: First Break, v. 26, no. 5, p. 109–113. Sudicky, E. A., and R. G. McLaren, 1992, The Laplace transform Galerkin technique for large-scale simulation of mass transport in discretely fractured porous formations: Water Resources Research, v. 28, no. 2, p. 499–514, doi:10.1029/91WR02560. Unsal, E., S. K. Matthäi, and M. J. Blunt, in press, Simulation of multiphase flow in fractured reservoirs using a fracture-only model with transfer functions: Computational Geosciences. van Duijn, C. J., and M. J. de Neef, 1998, Similarity solution for capillary redistribution of two phases in a porous medium with a single discontinuity: Advances in Water Resources, v. 21, p. 451–461. Wu, H., and D. D. Pollard, 1995, An experimental study of the relationship between joint spacing and layer thickness: Journal of Structural Geology, v. 17, p. 887–905, doi:10.1016/0191-8141 (94)00099-L. Zhang, X., and D. J. Sanderson, 1996, Effects of stress on the twodimensional permeability tensor of natural fracture networks: Geophysical Journal International, v. 125, p. 912–934, doi:10 .1111/j.1365-246X.1996.tb06034.x. Zimmermann, R. W., T. Hadgu, and G. S. Bodvarsson, 1996, A new lumped-parameter model for flow in unsaturated dual-porosity media: Advances in Water Resources, v. 19, no. 5, p. 317–327.