# Comparison of Flow Solutions for Naturally Fractured Reservoirs Using Complex Analysis Methods (CAM) and Embedded Discrete Fracture Models (EDFM): Fundamental Design Differences and Improved Scaling Method

Authors: Aaditya Khanal; Ruud Weijermars
Journal: Geofluids
 (2020)

---

## Abstract

The present study compares flow paths in reservoirs with natural fractures, solved with Complex Analysis Methods (CAM), to those solved with Embedded Discrete Fracture Models (EDFM). One aim is to define scaling rules for the strength (flux) of the discrete natural fractures used in CAM models, which was previously theoretically defined based on the expected flow distortion. A major hurdle for quantitative benchmarks of CAM with EDFM results is that each of the two methods accounts for natural fractures with different assumptions and input parameters. For example, EDFM scales the permeability of the natural fractures based on a cubic equation, while CAM uses a flux strength. The results from CAM and EDFM are used to scale the flux strength of the natural fractures and improve the equivalent permeability contrast estimation for CAM. The permeability contrast for CAM is calculated from the ratio of the enhanced velocity inside natural fractures to the unperturbed matrix fluid velocity. A significant advantage of flow and pressure models based on CAM is the high resolution without complex gridding. Particle tracking results are presented for fractures with different hydraulic conductivity ranging from highly permeable to impervious.


---

## Body

## 1. Introduction

The flow of fluids in naturally fractured reservoirs is highly influenced by the permeability, porosity, density, orientation, and several other features of discrete natural fractures. When such natural fractures have an enhanced permeability, they become highly conductive. The highly conductive natural fractures may alter the flow path of fluids by altering the local state of pressure and flow rates in the reservoir. The permeability contrast of natural fractures with the matrix also changes the shape of the drained rock volume by providing preferential flow paths to the trapped oil and gas fluids further away from a well [1, 2]. Natural fractures may also be reactivated during hydraulic fracturing, as predicted by fracture propagation models [3, 4] and microseismic events [5], thereby distorting the flow path. Unlike hydraulic fractures, natural fracture networks do not directly drain the reservoir as there is no pressure sink [6]. Nonetheless, prior studies have also shown that natural fractures may enhance pressure communication and flow interference between adjoining wells [1, 5].In the present study, we compare the results of CAM with EDFM for natural fractures oriented at different directions with respect to the principal direction of fluid flow. The three-term CAM algorithm [7] developed to trace the particle path deflection by natural fractures oriented at any angle with respect to the far-field flow is used. The outline of this paper is as follows. Section 2 provides a brief review of flow modeling tools for fractured porous media. Section 3 presents a background on the flux strength variable used in CAM, which relates to the hydraulic conductivity of natural fractures (Section 3.1). Section 3 also shows the scaling of the strength variable by comparing the pressure contour plots for natural fractures with different flux strengths, using CAM (Section 3.2), with EDFM results (Section 3.3). Section 4 shows CAM results (pressure contours and particle path solutions) for flow in naturally fractured reservoirs with arbitrary single fractures (Section 4.1). The CAM code can equally well account for impervious fractures that do not attract but impede fluid flow (Section 4.2). Section 5 highlights the fast computation time of CAM models (Section 5.1) and presents the results from a first attempt to compare CAM and EDFM models for flow across a reservoir section with multiple randomly oriented fractures (Section 5.2). A discussion follows in Section 6, and conclusions are given in Section 7.


## 2. Review of Flow Models for Naturally Fractured Reservoirs

Accurate simulation of fluid flow in naturally fractured porous media is a major challenge in subsurface reservoir engineering. The dimensions of the natural fractures are several orders of magnitude smaller than the dimensions of the reservoir, which adds to the complexity of accurately representing the natural fractures. We briefly summarize some of the approaches for modeling natural fractures below.


### 2.1. Prior Approaches

#### 2.1.1. Dual Continuum Models

Naturally fractured reservoirs were first modeled using a dual continuum approach formulation [8, 9]. In a dual-porosity model, natural fractures are represented by homogeneous and isotropic matrix blocks separated by orthogonal uniform natural fractures. The dual-porosity model assumes that the flow of fluid stored in a noncommunicating matrix occurs through the fractures. The dual-porosity model was later modified [10, 11] to account for other complex behavior seen in naturally fractured reservoirs. Dual-porosity models are still used today to model naturally fractured reservoirs. The merits are its relative simplicity and computational efficiency as compared to other discrete fracture and fracture network models.However, the dual-porosity model is not accurate for cases where the fracture geometry is complex and asymmetric, as is the case for hydraulically and naturally fractured unconventional reservoirs [12]. Modifications such as multiple interacting continua [13], time-dependent shape factors [14], and explicit parameterization of time-dependent transfer functions [15, 16] have been proposed to tackle the shortcomings of the multiporosity model. Despite these modifications, multiporosity models cannot explicitly account for the density and orientations of the natural fractures, which leads to unrealistic results. The shape factor and transfer function may not fully capture the complex flow behavior due to detailed pressure and fluid saturation gradients in naturally fractured reservoirs [17]. Multicontinuum models do not make any explicit geometric distinction among matrix, fractures, and fracture intersection. Implicit representation of fractures and matrix in such models needs the flow to be represented by upscaled quantities [12].In contrast to multiporosity/multicontinuum models, the explicit numerical modeling of fractures in discrete fracture models is computationally intensive but conceptually simpler than the implicit (multicontinuum) models. The discrete representation of fractures can be broadly categorized into four principal groups, (1) Discrete Fracture Network (DFN), (2) Discrete Fracture-Matrix (DFM) model, (3) EDFM, and (4) other gridded solution methods [12]. They are briefly reviewed below.


#### 2.1.2. DFN

In a DFN model, the matrix is assumed to be impermeable, and the flow is expected to occur only through the discrete fracture networks. DFN models consider fluid flow and transport processes in a fractured rock through a system of connected natural fractures. The DFN method is useful for studying fluid flow and mass transport in the fractured rocks for which an equivalent continuum model is difficult to establish or not applicable. It can also be used to derive the equivalent continuum flow and transport properties in the fractured rock for subsequent use in faster, upscaled (but implicit) reservoir models [18, 19]. In a DFN model, the storage and flow of fluids occur only through the fracture networks, which is suitable for modeling low-permeability and low-porosity fractured media. For a low-porosity/low-permeability system with many dominant natural fractures, the continuum approximation may not be entirely valid as the flow through the matrix is assumed to be negligible compared to the fractures. The DFN models may also be used to perform large-scale simulations where the fractured reservoir properties need to be approximated through upscaling and homogenization into equivalent permeability tensors [20].


#### 2.1.3. DFM

In a DFM model, the fractures are modeled as lower-dimensional interfaces embedded in the rock matrix. The DFM model reduces the loss of accuracy due to upscaling by introducing realistic geometrical complexities. In a DFM model, the fluid resides in both porous matrix and explicit fractures, but the smaller fractures are integrated into the matrix with appropriate upscaling. The DFM model is suitable for reservoirs with several natural fractures where only a few dominant fractures contribute to fluid storage and flow. The upscaling of matrix permeability to account for nondominant fractures reduces the complexity and computational cost during meshing without foregoing the accuracy. The selection of nondominant fractures integrated into the matrix is usually based on the dimensions of the fractures [21].


#### 2.1.4. EDFM

EDFM uses nonconforming grids with respect to fracture-matrix connections (introduced by [22]) and is an extension of the classical DFM model. EDFM uses a hybrid approach, where the dual-porosity model is used for the smaller and medium fractures, and DFN is used to model the larger fractures [22]. Flows within the matrix and the fractures are proportionated by the pressure difference between them and are discretized separately [23]. EDFM allows for complex fractures to be implemented in conventionally structured matrix grids without the need for local grid refinement (LGR) in the vicinity of the fractures [19]. Other models, such as projection-based EDFM (pEDFM), have also been proposed to improve the traditional EDFM [24, 25].


#### 2.1.5. Other Gridded Methods

Other methods such as the extended finite-element method or XFEM [26, 27] and conforming mesh using triangles and Voronoi grids [28, 29] can also be used to discretize naturally fractured reservoirs. These advanced discretization techniques capture the discontinuity of pressures across the fracture surface while preserving acceptable resolution of the near-fracture dynamics [26].All of the concurrent numerical methods (multicontinuum, DFN, DFM, EDFM, and XFEM) use discretization or meshing as a pivotal step to simulate the flow of fluids through the naturally fractured reservoirs. The discretization may require refinements of the gridding to account for the heterogeneities, interaction of the fracture/matrix system, and flow within the fractures/matrix [30]. Multicontinuum models are discretized by using finite difference where several values for physical parameters are assigned to each medium. For discrete fracture models, finite-element methods are primarily utilized to model the discrete fractures. However, a significant drawback of such advanced discretization schemes is the computational complexity and difficulty in accurately representing the prototype with a finite number of grid blocks [12, 26, 31].Efficient meshing/gridding is the biggest bottleneck to reduce the computation time of discrete numerical methods due to the inherent geometric complexity of fracture networks [12]. In addition, some models, such as EDFM, are only valid for high-permeability fractures and cannot model impervious or low-permeability fractures (for example, due to cementation or clay decay in fracture zones) [25]. Recently, analytically solvable models using Green’s function for gas flow in complex fracture networks have been presented by Marder et al. [32], which was also numerically tested with Barnett shale reservoir properties by Eftekhari et al. [33]. The existing numerical and analytical models are powerful tools with several strengths and weaknesses.


### 2.2. New Approach with Semianalytical CAM Models

In this study, we present an alternative method (CAM), with a low computational load that can accurately model and visualize the flow in various kinds of naturally fractured reservoirs. Traditionally, CAM uses potential and stream functions to describe the physical transport of particles in basic flow fields [34]. The basic flows can be combined or superposed with each other to describe complex flows that occur in groundwater and oil reservoirs. Increasingly complex flows can be combined such that solutions satisfy the Cauchy-Riemann differential equation. The functions (complex potentials) which satisfy the Cauchy-Riemann differential equation define the two-dimensional flow of incompressible and irrotational fluids. The basic algorithms for potential theory have been extensively described in our earlier publications [4, 34, 35] and fluid-mechanics literature [35–38]. Models based on CAM have been previously used to model fluid flow in hydraulically fractured reservoirs [39–41]. This paper applies various complex potentials based on areal doublets to model the flow of fluids in natural fractures [7, 17].A significant strength of models based on CAM is the ability to solve for the flow equations without any gridding. Discrete fracture models rely on the application of unstructured grids, which increases the computational complexity and makes the real-field applications challenging [25]. CAM algorithms do not require extensive gridding or meshing, which enables the modeling of heterogeneous reservoirs with numerous discrete fractures. Consequently, CAM algorithms are computationally efficient and offer high resolution, which is especially beneficial for modeling flow in unconventional oil and gas reservoirs that involve many hydraulic and natural fractures. CAM models of flow in reservoirs involving multiple wells, hydraulic fractures acting as pressure sinks, and impermeable fractures and faults have been validated as producing accurate results [42, 43]. CAM particle paths closely matched with those obtained by independent methods (e.g., Eclipse) [42, 43].CAM models applied to naturally fractured reservoirs have been presented elsewhere [1, 2, 17, 41], but no benchmarks of results against other methods have been presented yet. The present study is aimed at filling (at least a part of) that gap. A previous study has already pointed out that hydraulic fractures (connected to a wellbore) act as pressure sinks and behave differently from natural fractures [6], assuming the natural fractures of concern are not connected to a wellbore or a hydraulic fracture. Connected natural fractures behave like an extension of the hydraulic fracture network.The key algorithm used to model natural fractures in CAM was first derived by superposing areal doublet solutions [17], which are accurate for flow through fractures aligned with a far-field flow and can model multiple fractures with different flux strengths. The natural fracture element [17] was recently augmented [7] to accommodate the particle paths for fractures oriented at a large angle to the far-field flow. As the algorithms based on CAM are multivalued in certain branch-cut locations [44], the augmented solution [7] also needed to circumvent the branch cuts to avoid discontinuity in pressure (potential function) plots. The solution was augmented by gradually superposing two areal doublets with transformed coordinates, based on the angle with the far-field flow. The rotated fracture element results in the correct particle paths, even when the fractures occur perpendicular to the principal flow direction. For the intermediate cases, when the angles between the direction of flow and the areal doublet range between 0° and 90°, the particle paths are solved by the superposition of the original element and the rotated element [7]. The key algorithms used in this study are summarized in Section 3 and Appendix A.


## 2.1. Prior Approaches

### 2.1.1. Dual Continuum Models

Naturally fractured reservoirs were first modeled using a dual continuum approach formulation [8, 9]. In a dual-porosity model, natural fractures are represented by homogeneous and isotropic matrix blocks separated by orthogonal uniform natural fractures. The dual-porosity model assumes that the flow of fluid stored in a noncommunicating matrix occurs through the fractures. The dual-porosity model was later modified [10, 11] to account for other complex behavior seen in naturally fractured reservoirs. Dual-porosity models are still used today to model naturally fractured reservoirs. The merits are its relative simplicity and computational efficiency as compared to other discrete fracture and fracture network models.However, the dual-porosity model is not accurate for cases where the fracture geometry is complex and asymmetric, as is the case for hydraulically and naturally fractured unconventional reservoirs [12]. Modifications such as multiple interacting continua [13], time-dependent shape factors [14], and explicit parameterization of time-dependent transfer functions [15, 16] have been proposed to tackle the shortcomings of the multiporosity model. Despite these modifications, multiporosity models cannot explicitly account for the density and orientations of the natural fractures, which leads to unrealistic results. The shape factor and transfer function may not fully capture the complex flow behavior due to detailed pressure and fluid saturation gradients in naturally fractured reservoirs [17]. Multicontinuum models do not make any explicit geometric distinction among matrix, fractures, and fracture intersection. Implicit representation of fractures and matrix in such models needs the flow to be represented by upscaled quantities [12].In contrast to multiporosity/multicontinuum models, the explicit numerical modeling of fractures in discrete fracture models is computationally intensive but conceptually simpler than the implicit (multicontinuum) models. The discrete representation of fractures can be broadly categorized into four principal groups, (1) Discrete Fracture Network (DFN), (2) Discrete Fracture-Matrix (DFM) model, (3) EDFM, and (4) other gridded solution methods [12]. They are briefly reviewed below.


### 2.1.2. DFN

In a DFN model, the matrix is assumed to be impermeable, and the flow is expected to occur only through the discrete fracture networks. DFN models consider fluid flow and transport processes in a fractured rock through a system of connected natural fractures. The DFN method is useful for studying fluid flow and mass transport in the fractured rocks for which an equivalent continuum model is difficult to establish or not applicable. It can also be used to derive the equivalent continuum flow and transport properties in the fractured rock for subsequent use in faster, upscaled (but implicit) reservoir models [18, 19]. In a DFN model, the storage and flow of fluids occur only through the fracture networks, which is suitable for modeling low-permeability and low-porosity fractured media. For a low-porosity/low-permeability system with many dominant natural fractures, the continuum approximation may not be entirely valid as the flow through the matrix is assumed to be negligible compared to the fractures. The DFN models may also be used to perform large-scale simulations where the fractured reservoir properties need to be approximated through upscaling and homogenization into equivalent permeability tensors [20].


### 2.1.3. DFM

In a DFM model, the fractures are modeled as lower-dimensional interfaces embedded in the rock matrix. The DFM model reduces the loss of accuracy due to upscaling by introducing realistic geometrical complexities. In a DFM model, the fluid resides in both porous matrix and explicit fractures, but the smaller fractures are integrated into the matrix with appropriate upscaling. The DFM model is suitable for reservoirs with several natural fractures where only a few dominant fractures contribute to fluid storage and flow. The upscaling of matrix permeability to account for nondominant fractures reduces the complexity and computational cost during meshing without foregoing the accuracy. The selection of nondominant fractures integrated into the matrix is usually based on the dimensions of the fractures [21].


### 2.1.4. EDFM

EDFM uses nonconforming grids with respect to fracture-matrix connections (introduced by [22]) and is an extension of the classical DFM model. EDFM uses a hybrid approach, where the dual-porosity model is used for the smaller and medium fractures, and DFN is used to model the larger fractures [22]. Flows within the matrix and the fractures are proportionated by the pressure difference between them and are discretized separately [23]. EDFM allows for complex fractures to be implemented in conventionally structured matrix grids without the need for local grid refinement (LGR) in the vicinity of the fractures [19]. Other models, such as projection-based EDFM (pEDFM), have also been proposed to improve the traditional EDFM [24, 25].


### 2.1.5. Other Gridded Methods

Other methods such as the extended finite-element method or XFEM [26, 27] and conforming mesh using triangles and Voronoi grids [28, 29] can also be used to discretize naturally fractured reservoirs. These advanced discretization techniques capture the discontinuity of pressures across the fracture surface while preserving acceptable resolution of the near-fracture dynamics [26].All of the concurrent numerical methods (multicontinuum, DFN, DFM, EDFM, and XFEM) use discretization or meshing as a pivotal step to simulate the flow of fluids through the naturally fractured reservoirs. The discretization may require refinements of the gridding to account for the heterogeneities, interaction of the fracture/matrix system, and flow within the fractures/matrix [30]. Multicontinuum models are discretized by using finite difference where several values for physical parameters are assigned to each medium. For discrete fracture models, finite-element methods are primarily utilized to model the discrete fractures. However, a significant drawback of such advanced discretization schemes is the computational complexity and difficulty in accurately representing the prototype with a finite number of grid blocks [12, 26, 31].Efficient meshing/gridding is the biggest bottleneck to reduce the computation time of discrete numerical methods due to the inherent geometric complexity of fracture networks [12]. In addition, some models, such as EDFM, are only valid for high-permeability fractures and cannot model impervious or low-permeability fractures (for example, due to cementation or clay decay in fracture zones) [25]. Recently, analytically solvable models using Green’s function for gas flow in complex fracture networks have been presented by Marder et al. [32], which was also numerically tested with Barnett shale reservoir properties by Eftekhari et al. [33]. The existing numerical and analytical models are powerful tools with several strengths and weaknesses.


## 2.1.1. Dual Continuum Models

Naturally fractured reservoirs were first modeled using a dual continuum approach formulation [8, 9]. In a dual-porosity model, natural fractures are represented by homogeneous and isotropic matrix blocks separated by orthogonal uniform natural fractures. The dual-porosity model assumes that the flow of fluid stored in a noncommunicating matrix occurs through the fractures. The dual-porosity model was later modified [10, 11] to account for other complex behavior seen in naturally fractured reservoirs. Dual-porosity models are still used today to model naturally fractured reservoirs. The merits are its relative simplicity and computational efficiency as compared to other discrete fracture and fracture network models.However, the dual-porosity model is not accurate for cases where the fracture geometry is complex and asymmetric, as is the case for hydraulically and naturally fractured unconventional reservoirs [12]. Modifications such as multiple interacting continua [13], time-dependent shape factors [14], and explicit parameterization of time-dependent transfer functions [15, 16] have been proposed to tackle the shortcomings of the multiporosity model. Despite these modifications, multiporosity models cannot explicitly account for the density and orientations of the natural fractures, which leads to unrealistic results. The shape factor and transfer function may not fully capture the complex flow behavior due to detailed pressure and fluid saturation gradients in naturally fractured reservoirs [17]. Multicontinuum models do not make any explicit geometric distinction among matrix, fractures, and fracture intersection. Implicit representation of fractures and matrix in such models needs the flow to be represented by upscaled quantities [12].In contrast to multiporosity/multicontinuum models, the explicit numerical modeling of fractures in discrete fracture models is computationally intensive but conceptually simpler than the implicit (multicontinuum) models. The discrete representation of fractures can be broadly categorized into four principal groups, (1) Discrete Fracture Network (DFN), (2) Discrete Fracture-Matrix (DFM) model, (3) EDFM, and (4) other gridded solution methods [12]. They are briefly reviewed below.


## 2.1.2. DFN

In a DFN model, the matrix is assumed to be impermeable, and the flow is expected to occur only through the discrete fracture networks. DFN models consider fluid flow and transport processes in a fractured rock through a system of connected natural fractures. The DFN method is useful for studying fluid flow and mass transport in the fractured rocks for which an equivalent continuum model is difficult to establish or not applicable. It can also be used to derive the equivalent continuum flow and transport properties in the fractured rock for subsequent use in faster, upscaled (but implicit) reservoir models [18, 19]. In a DFN model, the storage and flow of fluids occur only through the fracture networks, which is suitable for modeling low-permeability and low-porosity fractured media. For a low-porosity/low-permeability system with many dominant natural fractures, the continuum approximation may not be entirely valid as the flow through the matrix is assumed to be negligible compared to the fractures. The DFN models may also be used to perform large-scale simulations where the fractured reservoir properties need to be approximated through upscaling and homogenization into equivalent permeability tensors [20].


## 2.1.3. DFM

In a DFM model, the fractures are modeled as lower-dimensional interfaces embedded in the rock matrix. The DFM model reduces the loss of accuracy due to upscaling by introducing realistic geometrical complexities. In a DFM model, the fluid resides in both porous matrix and explicit fractures, but the smaller fractures are integrated into the matrix with appropriate upscaling. The DFM model is suitable for reservoirs with several natural fractures where only a few dominant fractures contribute to fluid storage and flow. The upscaling of matrix permeability to account for nondominant fractures reduces the complexity and computational cost during meshing without foregoing the accuracy. The selection of nondominant fractures integrated into the matrix is usually based on the dimensions of the fractures [21].


## 2.1.4. EDFM

EDFM uses nonconforming grids with respect to fracture-matrix connections (introduced by [22]) and is an extension of the classical DFM model. EDFM uses a hybrid approach, where the dual-porosity model is used for the smaller and medium fractures, and DFN is used to model the larger fractures [22]. Flows within the matrix and the fractures are proportionated by the pressure difference between them and are discretized separately [23]. EDFM allows for complex fractures to be implemented in conventionally structured matrix grids without the need for local grid refinement (LGR) in the vicinity of the fractures [19]. Other models, such as projection-based EDFM (pEDFM), have also been proposed to improve the traditional EDFM [24, 25].


## 2.1.5. Other Gridded Methods

Other methods such as the extended finite-element method or XFEM [26, 27] and conforming mesh using triangles and Voronoi grids [28, 29] can also be used to discretize naturally fractured reservoirs. These advanced discretization techniques capture the discontinuity of pressures across the fracture surface while preserving acceptable resolution of the near-fracture dynamics [26].All of the concurrent numerical methods (multicontinuum, DFN, DFM, EDFM, and XFEM) use discretization or meshing as a pivotal step to simulate the flow of fluids through the naturally fractured reservoirs. The discretization may require refinements of the gridding to account for the heterogeneities, interaction of the fracture/matrix system, and flow within the fractures/matrix [30]. Multicontinuum models are discretized by using finite difference where several values for physical parameters are assigned to each medium. For discrete fracture models, finite-element methods are primarily utilized to model the discrete fractures. However, a significant drawback of such advanced discretization schemes is the computational complexity and difficulty in accurately representing the prototype with a finite number of grid blocks [12, 26, 31].Efficient meshing/gridding is the biggest bottleneck to reduce the computation time of discrete numerical methods due to the inherent geometric complexity of fracture networks [12]. In addition, some models, such as EDFM, are only valid for high-permeability fractures and cannot model impervious or low-permeability fractures (for example, due to cementation or clay decay in fracture zones) [25]. Recently, analytically solvable models using Green’s function for gas flow in complex fracture networks have been presented by Marder et al. [32], which was also numerically tested with Barnett shale reservoir properties by Eftekhari et al. [33]. The existing numerical and analytical models are powerful tools with several strengths and weaknesses.


## 2.2. New Approach with Semianalytical CAM Models

In this study, we present an alternative method (CAM), with a low computational load that can accurately model and visualize the flow in various kinds of naturally fractured reservoirs. Traditionally, CAM uses potential and stream functions to describe the physical transport of particles in basic flow fields [34]. The basic flows can be combined or superposed with each other to describe complex flows that occur in groundwater and oil reservoirs. Increasingly complex flows can be combined such that solutions satisfy the Cauchy-Riemann differential equation. The functions (complex potentials) which satisfy the Cauchy-Riemann differential equation define the two-dimensional flow of incompressible and irrotational fluids. The basic algorithms for potential theory have been extensively described in our earlier publications [4, 34, 35] and fluid-mechanics literature [35–38]. Models based on CAM have been previously used to model fluid flow in hydraulically fractured reservoirs [39–41]. This paper applies various complex potentials based on areal doublets to model the flow of fluids in natural fractures [7, 17].A significant strength of models based on CAM is the ability to solve for the flow equations without any gridding. Discrete fracture models rely on the application of unstructured grids, which increases the computational complexity and makes the real-field applications challenging [25]. CAM algorithms do not require extensive gridding or meshing, which enables the modeling of heterogeneous reservoirs with numerous discrete fractures. Consequently, CAM algorithms are computationally efficient and offer high resolution, which is especially beneficial for modeling flow in unconventional oil and gas reservoirs that involve many hydraulic and natural fractures. CAM models of flow in reservoirs involving multiple wells, hydraulic fractures acting as pressure sinks, and impermeable fractures and faults have been validated as producing accurate results [42, 43]. CAM particle paths closely matched with those obtained by independent methods (e.g., Eclipse) [42, 43].CAM models applied to naturally fractured reservoirs have been presented elsewhere [1, 2, 17, 41], but no benchmarks of results against other methods have been presented yet. The present study is aimed at filling (at least a part of) that gap. A previous study has already pointed out that hydraulic fractures (connected to a wellbore) act as pressure sinks and behave differently from natural fractures [6], assuming the natural fractures of concern are not connected to a wellbore or a hydraulic fracture. Connected natural fractures behave like an extension of the hydraulic fracture network.The key algorithm used to model natural fractures in CAM was first derived by superposing areal doublet solutions [17], which are accurate for flow through fractures aligned with a far-field flow and can model multiple fractures with different flux strengths. The natural fracture element [17] was recently augmented [7] to accommodate the particle paths for fractures oriented at a large angle to the far-field flow. As the algorithms based on CAM are multivalued in certain branch-cut locations [44], the augmented solution [7] also needed to circumvent the branch cuts to avoid discontinuity in pressure (potential function) plots. The solution was augmented by gradually superposing two areal doublets with transformed coordinates, based on the angle with the far-field flow. The rotated fracture element results in the correct particle paths, even when the fractures occur perpendicular to the principal flow direction. For the intermediate cases, when the angles between the direction of flow and the areal doublet range between 0° and 90°, the particle paths are solved by the superposition of the original element and the rotated element [7]. The key algorithms used in this study are summarized in Section 3 and Appendix A.


## 3. Scaling of CAM Fracture Strength to Permeability with EDFM

This study presents a comparison and scaling rules for models based on complex analysis methods (CAM) with EDFM models. For CAM-based algorithms, the flow intensity is scaled by the strength variable. The units of the strength depend on the complex potential of the flow element. For example, the strength of a vertical well has units of m2/s, the strength of a hydraulic fracture has units of m3/s, and the strength of an areal doublet/natural fracture element has units of m4/s. The strength variable can be positive (e.g., injectors) or negative (e.g., producers). This section discusses the scaling of the strength variable, which is a key input parameter for models based on CAM.The permeability of a reservoir is one of the most important variables that determine the productivity and deliverability a well. Fundamentally, the permeability of a porous medium is the proportionality constant in Darcy’s law, which defines the relationship between the pressure gradient and the fluid flux (flow rate per unit area). Reservoir permeability is an intrinsic property of the porous medium. Although for any particular rock type a higher permeability generally correlates with higher porosity [45], these parameters affect the time of flight (TOF) of migrating fluids in opposite directions [42]. The time of flight decreases when the permeability increases and slows down when the permeability decreases. For permeability, the opposite occurs: flow speeds up when the porosity decreases (thus shortens TOF), and flow slows down when the porosity increases (thus lengthens TOF).Permeability for a reservoir can be estimated from well logs, using empirical models such as the Carman-Kozeny equation [45]. Advanced logging tools such as nuclear magnetic resonance (NMR) are also used to calculate the formation permeability applying the Timur-Coates model [46] and the Schlumberger-Doll-Research model [47], especially in reservoirs where the Carman-Kozeny model does not work well. Where production data and limited reservoir characteristics are available, history matching can estimate the reservoir permeability [40]. For a naturally fractured reservoir, permeability can be broadly divided into matrix and fracture permeability, which are both measured in Darcy. Unconventional reservoirs are characterized by heterogeneous geology where each feature, including the fracture and matrix permeability, is significantly different from one region to another. As the characterization of each of the fractures is difficult, flow is simulated by using the upscaled permeability for single and multicontinuum models. For discrete models, the permeability for each individual fracture is assigned based on applicable statistical distributions [48].


### 3.1. CAM Strength Scaling of Natural Fractures

From Darcy’s law, time-dependent flow rate,Vf, across a natural fracture of length L, due to the time-dependent pressure gradient ΔPft is defined as follows [17]:
(1)Vft=−kfμΔPftL,m⋅s−1,where kf/μ is the ratio of fracture permeability to fluid viscosity. Similarly, the time-dependent flow rate, Vm, across a section of the matrix with the same length due to the time-dependent pressure gradient ΔPmt is defined as follows [17]:
(2)Vmt=−kmμΔPmtL,m⋅s−1.Assuming that the pressure gradient across the natural fracture and adjoining matrix is equal (i.e.,ΔPmt=ΔPft), from equations (1) and (2), the permeability ratio kf:km=Rk can be calculated as follows [17]:
(3)VfVm=kfkm=Rk.According to equation (3), the fluid velocity can be used to scale the fracture permeability based on the known matrix permeability, and vice versa. When both fracture and matrix permeability are unknown, the permeability contrast Rk of the matrix and the natural fracture can be calculated from the ratio of the respective fluid velocities.In our previous study, the permeability ratioRk was linked to the primary input parameters of CAM [6], which are modified here to account for the superposition of Vf being the result of Vm and the superposed flux. Our prior study [6] computed the fracture strength (υf) as υf=RkVmHfWfLf, where Vm is the far-field velocity in the matrix, and Hf, Wf, and Lf are the height, width, and length of the natural fracture element, respectively. However, the velocity in the fracture, Vf, is due to a preexisting Vm plus an additional velocity component Vflux, superposed due to the fracture flux:
(4)Vf=Vm+Vflux.Substituting now the modified fracture strengthυf′=VfluxHfWfLf into equation (4) and using Rk yields the following:
(5)Rk−1=υf′VmHfWfLf.The strength of the natural fracture and the far-field velocity in the matrix are denoted byυf′ and Vm, respectively. Thus, if the defined permeability contrast is known, the strength of the natural fracture element to be used for a discrete natural fracture in CAM can be calculated as follows:
(6)υf′=Rk−1VmHfWfLf.In the following section, the permeability contrasts for CAM fractures are calculated from equation (3) by taking the ratio of the maximum velocity due to the natural fracture and the original far-field velocity. However, as seen later in Section 3.3, the permeability contrast (Rk) calculated from equation (3) underestimates the permeability contrast calculated from the traditional models, and a correction factor is required. A correction factor, ξ, is introduced in Section 3.3 (see equation (7)) to scale a required permeability contrast (kf/km) with a corresponding natural fracture strength (υf′) for CAM, using the pressure contours generated from EDFM to calibrate the correction factor empirically. In what follows, we simply use υf, in lieu of υf′, by dropping the apostrophe.


### 3.2. Determination of Permeability Contrast with Matrix and Fracture Using CAM

In this section, a reservoir space is considered with an arbitrary, uniform far-field flow velocity of1.117×10−7m/s (3.5 m/year) from left to right. Figure 1 shows the pressure (Pa) field for the reservoir due to the far-field flow, where the pressure contours are perpendicular to the flow direction of the fluid. The reservoir is assumed to have permeability and porosity of 10 mD and 10%, respectively. The porosity scales the far-field flow rate up to a net effective velocity of 1.117×10−6m/s. In the remainder of this study, the effective strength and the effective velocity are reported, accounting for the porosity of the reservoir.Figure 1
Pressure contour (Pa) for a reservoir with far-field flow with a net effective velocity of 35 m/year from left to right. Initial reservoir pressure is zero. Scale of pressure change is relative to the right-end boundary of frame viewed.Next, consider a reservoir model (Figure2) with a solitary natural fracture located centrally in the flow space. The reservoir and fluid properties for the naturally fractured reservoir are summarized in Table 1.Figure 2
(a) Pressure contours (in Pa) for a reservoir with far-field flow and a single natural fracture. The effective far-field flow velocity is1.12×10−6m/s from left to right, and effective fracture strength is 1.12×10−6m4/s. (b) Velocity magnitude for the reservoir in (a). The maximum velocity of 6.7×10−6m/s is observed at the center of the fracture. (c) The portion of the natural fracture zoomed in (marked by a square in (b)) to examine the velocity around the natural fracture. The effect of the natural fracture on velocity changes is limited to areas extremely close to the natural fracture. (d) The corresponding permeability contrast (Rk) along the y-axis at x=4.5, calculated using equation (3).
(a)(b)(c)(d)Table 1
Attributes for a model reservoir simulated with CAM.
Natural fracture attributesSymbolValueNatural fracture width (m)W0.04Natural fracture length (m)L5Natural fracture height (m)H1Natural fracture angle to far-field flowα0Porosityn0.1Effective far-field flow rate (m/s)Vx1.117×10−6Effective natural fracture strength (m4/s)υf1.117×10−6Matrix permeability (mD)k10Viscosity (cP)μ1Figures2(a) and 2(b) show the pressure contours and the velocity field for the naturally fractured reservoir using the conditions of Table 1. The pressure contours (Figure 2(a)) are only mildly perturbed (compared to Figure 1), mainly near the natural fracture tips. Figure 2(b) shows that the velocity increases locally inside the natural fracture and slightly around the tips of the natural fracture. A small section of the natural fracture, marked by the square box in Figure 2(b), is maximized to closely examine the velocity in and near the fracture. The maximized portion (Figure 2(c)) shows that the maximum velocity at the center of the natural fracture is almost six times the original effective far-field velocity. Figure 2(c) also indicates that the increase in velocity outside of the natural fracture is negligible. The permeability contrast (Rk) between the matrix and the fracture is calculated to be 5.97 (equation (3)).The localized velocity changes across a natural fracture are further highlighted in Figure2(d), where the implied permeability contrast across the y-axis at x=4.5 is plotted. The cross-section (Figure 2(d)) shows that the natural fracture increases the velocity only within the natural fracture itself, with a negligible impact on the matrix velocity. The presence of a single natural fracture will have a negligible impact on the average or upscaled equivalent permeability of the reservoir. However, if the natural fracture density is high due to the presence of numerous natural fractures, the upscaled permeability for a reservoir may significantly increase (see upscaling in [6, 48]).This fluid velocity increase inside the natural fractures (Figure2(d)) is the primary reason for preferential flow paths and flow channeling due to fracture networks [49–52]. Multiple tracer transport studies on core samples of different length scales have shown that flow inside a fractured reservoir is highly heterogeneous [52–55]. Natural fractures may also result in fracture/well communication; thus, they need to be adequately accounted for while designing the infill wells and hydraulic fractures [1, 2]. The presence of natural fractures may alter the flow paths shifting the drained rock volume due to the local increase in fluid velocity (as shown in Figure 2(d)). For completeness, we refer to an earlier study, where the flow inside natural fractures was studied with a higher resolution than used in Figure 2(d). Although CAM is gridless, the plotting procedure is grid based and may falsely suggest a triangular-shaped flow profile in a narrow fracture if the solution grid chosen is overly coarse (for improved computational speed). When solved with sufficiently tight grid spacing, CAM-based velocity profiles inside natural fractures will be U-shaped [6].Figure3 varies the effective natural fracture strength to show the effect on pressure contours, using multipliers of 10, 20, and 30, as summarized in Table 2. Figure 3 shows that the increase in the effective strength of natural fractures is reflected in the enhanced curvature of the pressure contours near the fracture tips. For each of the sensitivity cases, the velocity profiles (not shown) resemble Figure 2(c). The maximum velocity occurs inside the natural fractures, and the velocity elsewhere in the reservoir is unaffected. Table 2 (second row) includes the maximum velocity contrast for each of the cases in Figure 3. Figures 3(b) and 3(c) have significantly higher effective natural fracture strengths compared to the base case of Figure 2, which results in the branch cuts around the fracture tips becoming more pronounced [44]. The increasing impact of branch cuts when the effective fracture strength increases is further illustrated in Appendix B.Figure 3
Pressure contours (in Pa) for a naturally fractured reservoir and a single natural fracture with a strength of (a)1.12×10−5, (b) 2.23×10−5, and (c) 3.35×10−5m4/s. The effective far-field velocity is 1.12×10−6m/s from left to right.
(a)(b)(c)Table 2
Inputs for the sensitivity of natural fracture strength.
QuantityFigure3(a)Figure3(b)Figure3(c)υf (m4/s)1.12×10−52.23×10−53.35×10−5Vf/Vm50.7100150


### 3.3. Scaling Natural Fracture Strength to Permeability Using Embedded Discrete Fracture Model (EDFM)

In this section, EDFM is used to scale the natural fracture strength to the units of permeability (mD). A simple fracture model based on the inputs in Table3 (and Figure 4) is used to generate the pressure contours for both the CAM and EDFM models. The results are visually inspected and iterated to generate closely matching pressure contours. Based on the comparison between the results from the two models, the strength for CAM is determined and scaled to generate the same pressure contours using EDFM for various permeability contrasts.Table 3
Natural fracture attributes for the simplified synthetic model (EDFM and CAM).
Natural fracture attributesSymbolValueNatural fracture width (m)W0.04Natural fracture length (m)L5Natural fracture height (m)H1Natural fracture angle to far-field flowα0°Boundary conditions (Pa)Left: 107; right: 0Fracture permeability (D)(a) 0.5; (b) 1; (c) 1.5Matrix permeability (D)k0.01Viscosity (cP)μ1Figure 4
Pressure contours (in Pa) generated from EDFM. The flow is driven by the Dirichlet boundary conditions of 107 Pa and 0 Pa at the right and left faces, respectively. The matrix permeability is 10 mD, and the permeability contrasts (kf/km) are (a) 50, (b) 100, and (c) 150.
(a)(b)(c)EDFM is a special form of a discrete fracture model (DFM) model (see Section2.1), introduced [22] to reduce the high computational cost associated with traditional DFM methods. EDFM defines fractures explicitly, as major fluid pathways, and benefits from independent definitions of the fracture and matrix grid. Thus, EDFM does not require a conforming mesh for the discrete fractures, which reduces the gridding complexity. Several authors have recently presented EDFM as a promising alternative to DFN and other upscaled single/multicontinuum models [25, 52, 53].Figure4 considers a single-phase flow in a 9×9m2 homogeneous domain where the matrix contains 225×225 grid cells, and one fracture occurs at the center of the domain. The natural fracture has a length and an aperture of 5 m and 0.04 m, respectively, and contains 50 grid cells with an average size of 0.1 m2. The flux for each (matrix-matrix and fracture-fracture) grid interface is defined by using the two-point-flux approximation (TPFA). The fracture cells are introduced into discrete systems through nonneighboring connections (NNC). The flow simulations are performed using MATLAB Reservoir Simulation Toolbox (MRST), a full-physics reservoir simulator, with EDFM or the hierarchical fracture model (HFM) [30]. The governing equations, formulations, and modeling techniques using EDFM in MRST are discussed in the literature [56–59].The flow in Figure4 is driven by the Dirichlet boundary conditions of 107 Pa and 0 Pa at the right and left faces of the flow domain, respectively. The matrix permeability is 0.01 D, and the fluid viscosity is 1 cP. The fracture permeability is set to the values of 0.5 D, 1 D, and 1.5 D representing the permeability contrast (Rk) of 50, 100, and 150, respectively. Table 3 summarizes the reservoir attributes, and Figure 4 shows the pressure contours for the flow domain for all the three cases.The pressure contours from EDFM (Figure4) show curving similar to the results from CAM (Figure 3). Several more pressure contour plots were generated using CAM, in addition to the results in Figure 3, to compare the results in Figure 4. Figures 5(a), 5(b), and 5(c) show the pressure contours from CAM for Vf/Vm values of 20, 40, and 60, respectively (shown in Figures 5(d), 5(e), and 5(f)).Figure 5
(a, b, and c) Pressure contours (in Pa) generated from CAM. The far-field flow moves from left to right with an effective velocity of1.12×10−6m/s. (d, e, and f) The respective ratio of fracture velocity (Vf) to matrix velocity (Vm) for fractures in (a), (b), and (c). The velocity contrast Vf/Vm is (a) 20, (b) 40, and (c) 60. The effective strength of the fractures for each case is (a) 1.12×10−5, (b) 2.23×10−5, and (c) 3.35×10−5m4/s, calculated from equation (7), where ξ=2.5.
(a)(b)(c)(d)(e)(f)The CAM pressure contours in Figures5(a), 5(b), and 5(c) show a better match with the pressure contours in Figure 4 (generated from EDFM). However, the permeability contrast for CAM (Figures 5(d), 5(e), and 5(f)), calculated from equation (3), which uses strength as the proxy for permeability, would be lower than the actual permeability contrast used in the EDFM model (Figure 4) by a factor of approximately 2.5. For instance, the Rk for the model generated from EDFM (Figure 4(a)) is 50, whereas the Rk for CAM is 20 (Figure 5(d)). Hence, an empirical correction factor ξ is introduced to scale the Rk calculated as the ratio of fracture and matrix permeability and to calculate the strength of the fractures υf, by using the modification of equation (6) as follows:
(7)υf=Rk−1ξVmHfWfLf,where Rk is the permeability contrast calculated from the ratio of fracture and matrix permeability in Darcy.Equation (7) facilitates flow modeling in naturally fractured porous media with CAM when the permeability contrast between the matrix and the fracture is known. The natural fracture strength, which is analogous to permeability, can be scaled using equation (7).


## 3.1. CAM Strength Scaling of Natural Fractures

From Darcy’s law, time-dependent flow rate,Vf, across a natural fracture of length L, due to the time-dependent pressure gradient ΔPft is defined as follows [17]:
(1)Vft=−kfμΔPftL,m⋅s−1,where kf/μ is the ratio of fracture permeability to fluid viscosity. Similarly, the time-dependent flow rate, Vm, across a section of the matrix with the same length due to the time-dependent pressure gradient ΔPmt is defined as follows [17]:
(2)Vmt=−kmμΔPmtL,m⋅s−1.Assuming that the pressure gradient across the natural fracture and adjoining matrix is equal (i.e.,ΔPmt=ΔPft), from equations (1) and (2), the permeability ratio kf:km=Rk can be calculated as follows [17]:
(3)VfVm=kfkm=Rk.According to equation (3), the fluid velocity can be used to scale the fracture permeability based on the known matrix permeability, and vice versa. When both fracture and matrix permeability are unknown, the permeability contrast Rk of the matrix and the natural fracture can be calculated from the ratio of the respective fluid velocities.In our previous study, the permeability ratioRk was linked to the primary input parameters of CAM [6], which are modified here to account for the superposition of Vf being the result of Vm and the superposed flux. Our prior study [6] computed the fracture strength (υf) as υf=RkVmHfWfLf, where Vm is the far-field velocity in the matrix, and Hf, Wf, and Lf are the height, width, and length of the natural fracture element, respectively. However, the velocity in the fracture, Vf, is due to a preexisting Vm plus an additional velocity component Vflux, superposed due to the fracture flux:
(4)Vf=Vm+Vflux.Substituting now the modified fracture strengthυf′=VfluxHfWfLf into equation (4) and using Rk yields the following:
(5)Rk−1=υf′VmHfWfLf.The strength of the natural fracture and the far-field velocity in the matrix are denoted byυf′ and Vm, respectively. Thus, if the defined permeability contrast is known, the strength of the natural fracture element to be used for a discrete natural fracture in CAM can be calculated as follows:
(6)υf′=Rk−1VmHfWfLf.In the following section, the permeability contrasts for CAM fractures are calculated from equation (3) by taking the ratio of the maximum velocity due to the natural fracture and the original far-field velocity. However, as seen later in Section 3.3, the permeability contrast (Rk) calculated from equation (3) underestimates the permeability contrast calculated from the traditional models, and a correction factor is required. A correction factor, ξ, is introduced in Section 3.3 (see equation (7)) to scale a required permeability contrast (kf/km) with a corresponding natural fracture strength (υf′) for CAM, using the pressure contours generated from EDFM to calibrate the correction factor empirically. In what follows, we simply use υf, in lieu of υf′, by dropping the apostrophe.


## 3.2. Determination of Permeability Contrast with Matrix and Fracture Using CAM

In this section, a reservoir space is considered with an arbitrary, uniform far-field flow velocity of1.117×10−7m/s (3.5 m/year) from left to right. Figure 1 shows the pressure (Pa) field for the reservoir due to the far-field flow, where the pressure contours are perpendicular to the flow direction of the fluid. The reservoir is assumed to have permeability and porosity of 10 mD and 10%, respectively. The porosity scales the far-field flow rate up to a net effective velocity of 1.117×10−6m/s. In the remainder of this study, the effective strength and the effective velocity are reported, accounting for the porosity of the reservoir.Figure 1
Pressure contour (Pa) for a reservoir with far-field flow with a net effective velocity of 35 m/year from left to right. Initial reservoir pressure is zero. Scale of pressure change is relative to the right-end boundary of frame viewed.Next, consider a reservoir model (Figure2) with a solitary natural fracture located centrally in the flow space. The reservoir and fluid properties for the naturally fractured reservoir are summarized in Table 1.Figure 2
(a) Pressure contours (in Pa) for a reservoir with far-field flow and a single natural fracture. The effective far-field flow velocity is1.12×10−6m/s from left to right, and effective fracture strength is 1.12×10−6m4/s. (b) Velocity magnitude for the reservoir in (a). The maximum velocity of 6.7×10−6m/s is observed at the center of the fracture. (c) The portion of the natural fracture zoomed in (marked by a square in (b)) to examine the velocity around the natural fracture. The effect of the natural fracture on velocity changes is limited to areas extremely close to the natural fracture. (d) The corresponding permeability contrast (Rk) along the y-axis at x=4.5, calculated using equation (3).
(a)(b)(c)(d)Table 1
Attributes for a model reservoir simulated with CAM.
Natural fracture attributesSymbolValueNatural fracture width (m)W0.04Natural fracture length (m)L5Natural fracture height (m)H1Natural fracture angle to far-field flowα0Porosityn0.1Effective far-field flow rate (m/s)Vx1.117×10−6Effective natural fracture strength (m4/s)υf1.117×10−6Matrix permeability (mD)k10Viscosity (cP)μ1Figures2(a) and 2(b) show the pressure contours and the velocity field for the naturally fractured reservoir using the conditions of Table 1. The pressure contours (Figure 2(a)) are only mildly perturbed (compared to Figure 1), mainly near the natural fracture tips. Figure 2(b) shows that the velocity increases locally inside the natural fracture and slightly around the tips of the natural fracture. A small section of the natural fracture, marked by the square box in Figure 2(b), is maximized to closely examine the velocity in and near the fracture. The maximized portion (Figure 2(c)) shows that the maximum velocity at the center of the natural fracture is almost six times the original effective far-field velocity. Figure 2(c) also indicates that the increase in velocity outside of the natural fracture is negligible. The permeability contrast (Rk) between the matrix and the fracture is calculated to be 5.97 (equation (3)).The localized velocity changes across a natural fracture are further highlighted in Figure2(d), where the implied permeability contrast across the y-axis at x=4.5 is plotted. The cross-section (Figure 2(d)) shows that the natural fracture increases the velocity only within the natural fracture itself, with a negligible impact on the matrix velocity. The presence of a single natural fracture will have a negligible impact on the average or upscaled equivalent permeability of the reservoir. However, if the natural fracture density is high due to the presence of numerous natural fractures, the upscaled permeability for a reservoir may significantly increase (see upscaling in [6, 48]).This fluid velocity increase inside the natural fractures (Figure2(d)) is the primary reason for preferential flow paths and flow channeling due to fracture networks [49–52]. Multiple tracer transport studies on core samples of different length scales have shown that flow inside a fractured reservoir is highly heterogeneous [52–55]. Natural fractures may also result in fracture/well communication; thus, they need to be adequately accounted for while designing the infill wells and hydraulic fractures [1, 2]. The presence of natural fractures may alter the flow paths shifting the drained rock volume due to the local increase in fluid velocity (as shown in Figure 2(d)). For completeness, we refer to an earlier study, where the flow inside natural fractures was studied with a higher resolution than used in Figure 2(d). Although CAM is gridless, the plotting procedure is grid based and may falsely suggest a triangular-shaped flow profile in a narrow fracture if the solution grid chosen is overly coarse (for improved computational speed). When solved with sufficiently tight grid spacing, CAM-based velocity profiles inside natural fractures will be U-shaped [6].Figure3 varies the effective natural fracture strength to show the effect on pressure contours, using multipliers of 10, 20, and 30, as summarized in Table 2. Figure 3 shows that the increase in the effective strength of natural fractures is reflected in the enhanced curvature of the pressure contours near the fracture tips. For each of the sensitivity cases, the velocity profiles (not shown) resemble Figure 2(c). The maximum velocity occurs inside the natural fractures, and the velocity elsewhere in the reservoir is unaffected. Table 2 (second row) includes the maximum velocity contrast for each of the cases in Figure 3. Figures 3(b) and 3(c) have significantly higher effective natural fracture strengths compared to the base case of Figure 2, which results in the branch cuts around the fracture tips becoming more pronounced [44]. The increasing impact of branch cuts when the effective fracture strength increases is further illustrated in Appendix B.Figure 3
Pressure contours (in Pa) for a naturally fractured reservoir and a single natural fracture with a strength of (a)1.12×10−5, (b) 2.23×10−5, and (c) 3.35×10−5m4/s. The effective far-field velocity is 1.12×10−6m/s from left to right.
(a)(b)(c)Table 2
Inputs for the sensitivity of natural fracture strength.
QuantityFigure3(a)Figure3(b)Figure3(c)υf (m4/s)1.12×10−52.23×10−53.35×10−5Vf/Vm50.7100150


## 3.3. Scaling Natural Fracture Strength to Permeability Using Embedded Discrete Fracture Model (EDFM)

In this section, EDFM is used to scale the natural fracture strength to the units of permeability (mD). A simple fracture model based on the inputs in Table3 (and Figure 4) is used to generate the pressure contours for both the CAM and EDFM models. The results are visually inspected and iterated to generate closely matching pressure contours. Based on the comparison between the results from the two models, the strength for CAM is determined and scaled to generate the same pressure contours using EDFM for various permeability contrasts.Table 3
Natural fracture attributes for the simplified synthetic model (EDFM and CAM).
Natural fracture attributesSymbolValueNatural fracture width (m)W0.04Natural fracture length (m)L5Natural fracture height (m)H1Natural fracture angle to far-field flowα0°Boundary conditions (Pa)Left: 107; right: 0Fracture permeability (D)(a) 0.5; (b) 1; (c) 1.5Matrix permeability (D)k0.01Viscosity (cP)μ1Figure 4
Pressure contours (in Pa) generated from EDFM. The flow is driven by the Dirichlet boundary conditions of 107 Pa and 0 Pa at the right and left faces, respectively. The matrix permeability is 10 mD, and the permeability contrasts (kf/km) are (a) 50, (b) 100, and (c) 150.
(a)(b)(c)EDFM is a special form of a discrete fracture model (DFM) model (see Section2.1), introduced [22] to reduce the high computational cost associated with traditional DFM methods. EDFM defines fractures explicitly, as major fluid pathways, and benefits from independent definitions of the fracture and matrix grid. Thus, EDFM does not require a conforming mesh for the discrete fractures, which reduces the gridding complexity. Several authors have recently presented EDFM as a promising alternative to DFN and other upscaled single/multicontinuum models [25, 52, 53].Figure4 considers a single-phase flow in a 9×9m2 homogeneous domain where the matrix contains 225×225 grid cells, and one fracture occurs at the center of the domain. The natural fracture has a length and an aperture of 5 m and 0.04 m, respectively, and contains 50 grid cells with an average size of 0.1 m2. The flux for each (matrix-matrix and fracture-fracture) grid interface is defined by using the two-point-flux approximation (TPFA). The fracture cells are introduced into discrete systems through nonneighboring connections (NNC). The flow simulations are performed using MATLAB Reservoir Simulation Toolbox (MRST), a full-physics reservoir simulator, with EDFM or the hierarchical fracture model (HFM) [30]. The governing equations, formulations, and modeling techniques using EDFM in MRST are discussed in the literature [56–59].The flow in Figure4 is driven by the Dirichlet boundary conditions of 107 Pa and 0 Pa at the right and left faces of the flow domain, respectively. The matrix permeability is 0.01 D, and the fluid viscosity is 1 cP. The fracture permeability is set to the values of 0.5 D, 1 D, and 1.5 D representing the permeability contrast (Rk) of 50, 100, and 150, respectively. Table 3 summarizes the reservoir attributes, and Figure 4 shows the pressure contours for the flow domain for all the three cases.The pressure contours from EDFM (Figure4) show curving similar to the results from CAM (Figure 3). Several more pressure contour plots were generated using CAM, in addition to the results in Figure 3, to compare the results in Figure 4. Figures 5(a), 5(b), and 5(c) show the pressure contours from CAM for Vf/Vm values of 20, 40, and 60, respectively (shown in Figures 5(d), 5(e), and 5(f)).Figure 5
(a, b, and c) Pressure contours (in Pa) generated from CAM. The far-field flow moves from left to right with an effective velocity of1.12×10−6m/s. (d, e, and f) The respective ratio of fracture velocity (Vf) to matrix velocity (Vm) for fractures in (a), (b), and (c). The velocity contrast Vf/Vm is (a) 20, (b) 40, and (c) 60. The effective strength of the fractures for each case is (a) 1.12×10−5, (b) 2.23×10−5, and (c) 3.35×10−5m4/s, calculated from equation (7), where ξ=2.5.
(a)(b)(c)(d)(e)(f)The CAM pressure contours in Figures5(a), 5(b), and 5(c) show a better match with the pressure contours in Figure 4 (generated from EDFM). However, the permeability contrast for CAM (Figures 5(d), 5(e), and 5(f)), calculated from equation (3), which uses strength as the proxy for permeability, would be lower than the actual permeability contrast used in the EDFM model (Figure 4) by a factor of approximately 2.5. For instance, the Rk for the model generated from EDFM (Figure 4(a)) is 50, whereas the Rk for CAM is 20 (Figure 5(d)). Hence, an empirical correction factor ξ is introduced to scale the Rk calculated as the ratio of fracture and matrix permeability and to calculate the strength of the fractures υf, by using the modification of equation (6) as follows:
(7)υf=Rk−1ξVmHfWfLf,where Rk is the permeability contrast calculated from the ratio of fracture and matrix permeability in Darcy.Equation (7) facilitates flow modeling in naturally fractured porous media with CAM when the permeability contrast between the matrix and the fracture is known. The natural fracture strength, which is analogous to permeability, can be scaled using equation (7).


## 4. Application of the Augmented Solution

An augmented CAM solution for the areal doublet was proposed to more accurately account for the refraction of particle paths across fractures not aligned with the far-field flow [7]. The augmented CAM solution was obtained by superposing two different complex potentials. The first complex potential which is superposed is the original areal doublet proposed by van Harmelen and Weijermars [17]. The second complex potential superposed to obtain the augmented solution is obtained by rotating the vertices of the first areal doublet [7]. The contribution of the two elements is tuned by the Sine function. Figures 6(a) and 6(b) show the particle paths (blue) and the time-of-flight contours (TOFCs, red) obtained from the original solution [17] and the augmented solution [7], respectively. For both cases, the far-field flow with an effective velocity of 3.12×10−8m/s flows from bottom to top. An areal doublet element with an effective strength of 3.12×10−6m4/s is assumed to be present at the center of the flow domain. The natural fracture is assumed to be 10 m long and 1 m wide. The fluid particles at the bottom of the flow domain are tracked for 30 years by using the Eulerian particle tracking method. The areal doublet from the original solution (Figure 6(a)) accelerates most of the fluid towards the left opening of the channel, even when the far-field flow is flowing from a different direction. However, Figure 6(b) shows that for the augmented solution, the fluid velocity is increased, but the refraction of the particle paths stays symmetrical to the far-field flow as expected.Figure 6
Particle path (blue) and TOFC (red) for (a) the original areal doublet solution [17] and (b) the augmented areal doublet solution [7]. Models in (a) and (b) both have a fracture strength of 3.12×10−6m4/s. Despite the high incidence angle, the particle paths for the augmented solution in (b) stay mostly aligned with the far-field flow, even inside the fracture zone (except near the fracture tips).
(a)(b)In the remainder of this section, the augmented solution [7] is used to generate pressure contours for flow channels (i.e., the natural fractures) with different fracture apertures and permeability contrast. Section 4.1 presents the pressure contours for highly conductive fractures with different aperture sizes. Section 4.2 presents the pressure contours for natural fractures with reduced permeability relative to the ambient matrix rock in the reservoir.


### 4.1. Pressure Contour Sensitivity to Fracture Aperture

Accurate representation of natural fracture dimensions and properties is essential to generate precise continuum and discrete fracture models [21, 60]. The fracture aperture or width is one such parameter that governs the fracture porosity and permeability and are constrained by fracture surface topography, shear displacement, and confining stress [61, 62]. The distribution of open fracture apertures in the subsurface is highly variable, which depends on fracture type, host lithology, degree of mineral fill/dissolution, and the in situ stress regime [63, 64]. Apertures are usually estimated by using different probability distribution functions, such as log-normal, power-law distribution function, and uniform distribution, due to the lack of available subsurface data [63–65]. In this section, we investigate the effect of fracture width and orientation on pressure contour patterns using a single natural fracture.Figure7 shows the pressure contours generated with the augmented CAM model for natural fractures with different apertures oriented at various angles. The fracture aperture for Case A (Figure 7(a)A, C, E) and Case B (Figure 7(b)B, D, F) is 0.04 m and 1 m, respectively. Other fracture attributes used for the CAM simulations of Figure 7 are summarized in Table 4.Figure 7
Pressure contours (in Pa) for fracture properties listed in Table4. Case A (a): fracture width is 0.04 m. Case B (b): fracture width increased to 1 m. The far-field flow with a net effective velocity of 3.12×10−7m/s flows from bottom to top. The effective strength of the natural fracture is 2.38×10−6m4/s. The natural fracture is oriented at angles of (A and B) 0°, (C and D) 45°, and (E and F) 90°, with respect to far-field flow direction.
(a)(b)Table 4
Fracture and reservoir attributes for fractures oriented at different angles.
Fracture/reservoir attributesSymbolValueNatural fracture width (m)WCaseA=0.04CaseB=1Natural fracture length (m)L10Natural fracture height (m)H1Natural fracture angle to far-field flowα90°Porosityn0.1Effective far-field flow rate (m/s)Vx3.12×10−7Effective natural fracture strength (m4/s)υf2.38×10−6Matrix permeability (mD)k10Viscosity (cP)μ1The pressure contour patterns for both Cases A and B are distorted near the fracture tips (Figures7(a)A, C, E and 7(b)B, D, F). The pressure contours for the fractures, which are not parallel to the direction of fluid flow (Figures 7(a)C, E and 7(b)D, F), show pressure jumps due to the integral effects of locally non-single-valued functions, which create branch cuts [44]. The branch cuts may have a significant effect on the pressure contours, especially towards the fracture tips, as shown by Figure 8 in Appendix B, where the fracture strength was increased further as compared to Figure 7(a).Figure 8
Pressure contours (in Pa) for fractures with the same properties as Case A (Figure7(a)) but with doubled effective fracture strength of 4.76×10−6m4/s.Figure9 shows the particle paths for Cases A and B corresponding to the model conditions of Figures 7(a) and 7(b). The particle paths represent the progressive movement of fluid over 30 years. The particle paths (blue) and the time-of-flight contours (TOFCs, red, spaced at 3 years) show that the fluid moves further in fractures with smaller apertures, which are otherwise identical to each other. Thus, a smaller aperture may promote flow channeling. However, if we were to use a scaling of the fracture strength according to equation (7) and keep Rk constant but adjust the flux strength υf in proportion to the fracture width Wf, then both Cases A and B would have the same time of flight (TOF). Henceforth, it can be misleading to simply state that a smaller fracture aperture promotes flow channeling. Flow channeling is foremost an effect due to the permeability contrast between the matrix and the fracture, as expressed in Rk. A larger Rk will promote fracture channeling and lead to a shorter TOFC. When Rk<1, the fracture becomes progressively impervious, leading to a longer TOFC for fluid traveling via the fracture. Also, when Rk=1, the fracture may physically exist, but its presence will not affect the flow paths.Figure 9
Particle paths (blue) and time-of-flight contours, spaced at 3 years (red) for cases shown in Figure7.
(a)(b)


### 4.2. Pressure Contours for Impervious Fractures

Natural fractures may either be highly conductive or poorly conductive relative to the matrix, depending on the mineralization of the pore structure in the fracture zone [66]. For example, even cemented or blocked fractures can still be critical to the fracture network and may promote preferential flow channeling [67]. The degree of cementation in natural fractures depends on the burial conditions; original fracture aperture; and the geochemical environment, reactivity, and composition of the fracture wall rocks [68]. Most natural fractures in shale formations such as the Barnett are observed to be filled with calcite or quartz cement [68]. The cemented natural fractures also interact with hydraulic fractures to impede and divert the fracture propagation path [69, 70].In this section, we use CAM to generate the pressure contours for cemented or blocking fractures, where we revisit the natural fracture in Case A (Figure7(a)) with a small fracture aperture of 0.04 m. All the reservoir and flow properties are kept constant, except for the strength of the natural fracture. The effective strength of the natural fracture is set to −9.6×10−7m4/s. The negative sign opposes the far-field flow to mimic the action of a blocking/cemented fracture. Figure 10 shows the pressure contours (Figures 10(a), 10(b), and 10(c)) and the particle paths (Figures 10(c), 10(d), and 10(e)) for a simple blocking—but still permeable—fracture. The pressure contours for such a permeable, blocking fracture (Figures 10(a), 10(b), and 10(c)) show the opposite behavior to the highly conductive fracture (Figure 7(a)A, B, C). In both cases, the pressure contours are distorted near the fracture tips. The particle paths (blue) and the TOFCs (red) in Figures 10(c), 10(d), and 10(e) are generated by tracking a limited number of fluid particles originally, at the bottom of the plot, for 30 years. Each TOFC shows the movement of the far-field flow after three years, for a total flow time of 30 years. The particle paths show that the blocking fracture (but still slightly pervious) slows down the fluid particles in its path. The TOFCs around the natural fracture are pulled back by approximately 2.5 m in each case. Figure 11 in Appendix B shows the effect of increasing the flow resistance by further reducing the strength of the fracture by a factor of 2.Figure 10
(a, b, and c) Pressure contours (in Pa) for cemented fractures with the same properties as Case A (Figure7(a), Table 4) except for effective fracture strength, which has a value of −9.6×10−7m4/s. The natural fractures are oriented at an angle of (a) 0°, (b) 45°, and (c) 90°, with respect to the far-field flow. (c, d, and e) Particle paths and TOFCs for the corresponding cases in (a), (b), and (c).
(a)(b)(c)(d)(e)(f)Figure 11
Pressure contours (in Pa) for cemented fractures with the same properties as Figure11 but with effective fracture strength reduced by a factor of 2 to a value of −1.90×10−6m4/s.


## 4.1. Pressure Contour Sensitivity to Fracture Aperture

Accurate representation of natural fracture dimensions and properties is essential to generate precise continuum and discrete fracture models [21, 60]. The fracture aperture or width is one such parameter that governs the fracture porosity and permeability and are constrained by fracture surface topography, shear displacement, and confining stress [61, 62]. The distribution of open fracture apertures in the subsurface is highly variable, which depends on fracture type, host lithology, degree of mineral fill/dissolution, and the in situ stress regime [63, 64]. Apertures are usually estimated by using different probability distribution functions, such as log-normal, power-law distribution function, and uniform distribution, due to the lack of available subsurface data [63–65]. In this section, we investigate the effect of fracture width and orientation on pressure contour patterns using a single natural fracture.Figure7 shows the pressure contours generated with the augmented CAM model for natural fractures with different apertures oriented at various angles. The fracture aperture for Case A (Figure 7(a)A, C, E) and Case B (Figure 7(b)B, D, F) is 0.04 m and 1 m, respectively. Other fracture attributes used for the CAM simulations of Figure 7 are summarized in Table 4.Figure 7
Pressure contours (in Pa) for fracture properties listed in Table4. Case A (a): fracture width is 0.04 m. Case B (b): fracture width increased to 1 m. The far-field flow with a net effective velocity of 3.12×10−7m/s flows from bottom to top. The effective strength of the natural fracture is 2.38×10−6m4/s. The natural fracture is oriented at angles of (A and B) 0°, (C and D) 45°, and (E and F) 90°, with respect to far-field flow direction.
(a)(b)Table 4
Fracture and reservoir attributes for fractures oriented at different angles.
Fracture/reservoir attributesSymbolValueNatural fracture width (m)WCaseA=0.04CaseB=1Natural fracture length (m)L10Natural fracture height (m)H1Natural fracture angle to far-field flowα90°Porosityn0.1Effective far-field flow rate (m/s)Vx3.12×10−7Effective natural fracture strength (m4/s)υf2.38×10−6Matrix permeability (mD)k10Viscosity (cP)μ1The pressure contour patterns for both Cases A and B are distorted near the fracture tips (Figures7(a)A, C, E and 7(b)B, D, F). The pressure contours for the fractures, which are not parallel to the direction of fluid flow (Figures 7(a)C, E and 7(b)D, F), show pressure jumps due to the integral effects of locally non-single-valued functions, which create branch cuts [44]. The branch cuts may have a significant effect on the pressure contours, especially towards the fracture tips, as shown by Figure 8 in Appendix B, where the fracture strength was increased further as compared to Figure 7(a).Figure 8
Pressure contours (in Pa) for fractures with the same properties as Case A (Figure7(a)) but with doubled effective fracture strength of 4.76×10−6m4/s.Figure9 shows the particle paths for Cases A and B corresponding to the model conditions of Figures 7(a) and 7(b). The particle paths represent the progressive movement of fluid over 30 years. The particle paths (blue) and the time-of-flight contours (TOFCs, red, spaced at 3 years) show that the fluid moves further in fractures with smaller apertures, which are otherwise identical to each other. Thus, a smaller aperture may promote flow channeling. However, if we were to use a scaling of the fracture strength according to equation (7) and keep Rk constant but adjust the flux strength υf in proportion to the fracture width Wf, then both Cases A and B would have the same time of flight (TOF). Henceforth, it can be misleading to simply state that a smaller fracture aperture promotes flow channeling. Flow channeling is foremost an effect due to the permeability contrast between the matrix and the fracture, as expressed in Rk. A larger Rk will promote fracture channeling and lead to a shorter TOFC. When Rk<1, the fracture becomes progressively impervious, leading to a longer TOFC for fluid traveling via the fracture. Also, when Rk=1, the fracture may physically exist, but its presence will not affect the flow paths.Figure 9
Particle paths (blue) and time-of-flight contours, spaced at 3 years (red) for cases shown in Figure7.
(a)(b)


## 4.2. Pressure Contours for Impervious Fractures

Natural fractures may either be highly conductive or poorly conductive relative to the matrix, depending on the mineralization of the pore structure in the fracture zone [66]. For example, even cemented or blocked fractures can still be critical to the fracture network and may promote preferential flow channeling [67]. The degree of cementation in natural fractures depends on the burial conditions; original fracture aperture; and the geochemical environment, reactivity, and composition of the fracture wall rocks [68]. Most natural fractures in shale formations such as the Barnett are observed to be filled with calcite or quartz cement [68]. The cemented natural fractures also interact with hydraulic fractures to impede and divert the fracture propagation path [69, 70].In this section, we use CAM to generate the pressure contours for cemented or blocking fractures, where we revisit the natural fracture in Case A (Figure7(a)) with a small fracture aperture of 0.04 m. All the reservoir and flow properties are kept constant, except for the strength of the natural fracture. The effective strength of the natural fracture is set to −9.6×10−7m4/s. The negative sign opposes the far-field flow to mimic the action of a blocking/cemented fracture. Figure 10 shows the pressure contours (Figures 10(a), 10(b), and 10(c)) and the particle paths (Figures 10(c), 10(d), and 10(e)) for a simple blocking—but still permeable—fracture. The pressure contours for such a permeable, blocking fracture (Figures 10(a), 10(b), and 10(c)) show the opposite behavior to the highly conductive fracture (Figure 7(a)A, B, C). In both cases, the pressure contours are distorted near the fracture tips. The particle paths (blue) and the TOFCs (red) in Figures 10(c), 10(d), and 10(e) are generated by tracking a limited number of fluid particles originally, at the bottom of the plot, for 30 years. Each TOFC shows the movement of the far-field flow after three years, for a total flow time of 30 years. The particle paths show that the blocking fracture (but still slightly pervious) slows down the fluid particles in its path. The TOFCs around the natural fracture are pulled back by approximately 2.5 m in each case. Figure 11 in Appendix B shows the effect of increasing the flow resistance by further reducing the strength of the fracture by a factor of 2.Figure 10
(a, b, and c) Pressure contours (in Pa) for cemented fractures with the same properties as Case A (Figure7(a), Table 4) except for effective fracture strength, which has a value of −9.6×10−7m4/s. The natural fractures are oriented at an angle of (a) 0°, (b) 45°, and (c) 90°, with respect to the far-field flow. (c, d, and e) Particle paths and TOFCs for the corresponding cases in (a), (b), and (c).
(a)(b)(c)(d)(e)(f)Figure 11
Pressure contours (in Pa) for cemented fractures with the same properties as Figure11 but with effective fracture strength reduced by a factor of 2 to a value of −1.90×10−6m4/s.


## 5. Comparison of CAM and EDFM Results

In the approach below, we attempt a first benchmark of CAM-based solutions for flow in naturally fractured porous media with EDFM. There are several differences in the model design parameters of CAM and EDFM that may impede a direct comparison of model results, as explained here. For example, for modeling purposes, a fractured porous medium may be represented by a primary permeability due to connected pores and a secondary permeability due to the fracture conduits [70]. In numerical models, the secondary permeability is commonly an open fracture with a pseudopermeability assigned, which is a value based on the fracture aperture using a cubic equation. Such open fractures are likely to have an assigned permeability, being several orders of magnitude larger than the primary permeability, which will have a major impact on the upscaled equivalent (or effective) permeability of a representative elementary volume [6]. In CAM, fractures are assigned a strength that scales the permeability contrast with the matrix. Despite these significant differences, we present a first attempt to benchmark the results from CAM as compared to EDFM (Section 5.1), followed by a demonstration of runtime results for CAM using different time steps (Section 5.2).


### 5.1. Model Description and Results

The model design and particle paths from the EDFM reference solutions of Shah et al. [56] were used as a starting point for comparison with CAM results (Figures 12(a) and 12(b)). The EDFM representation originated from a dual continuum model based on a square of unit dimensions transected by fractures with unstated apertures, and geometry as portrayed in Figure 12(a). The fracture-matrix permeability ratio Rk=kf/km which is constant for all fractures is 104. Figure 12(b) shows the flow paths generated from EDFM based on the boundary conditions in [56]. The boundary conditions in the benchmark of Figure 12(b) are a combination of two Neumann’s conditions (namely no-flow boundaries, such that the directional derivative vanishes) and two Dirichlet’s conditions (with directional derivatives normal to the left and right boundary being due to uniform pressures). Essentially, no-flow boundaries are imposed at the upper and lower boundaries of the 2D flow area studied, and constant, uniform pressures are maintained at the left and right boundaries.In our CAM model, a no-flow upper boundary and lower boundary can be simulated by a fracture channel, but this would require Schwarz-Christoffel’s and Schottky-Kleine’s prime function boundary mapping as was used in our models of bounded reservoirs [71]. However, CAM is most user-friendly when applied to unbound flow domains, which is why we assume a uniform flux entering the flow space of the unit square from the left boundary (akin to a uniform pressure) (Figure 12(b)). The fracture lengths and orientations are extracted from Figure 12(a). The width of the fractures is not given [56] and is assumed to be 0.2 mm (0.0002 m). Although the available field data for the width (or aperture) of a natural fracture are limited from subsurface observations, values of 0.01 to 10 mm are reported in the literature [64]. The permeability ratios are determined by scaling the fracture flow strengths following the procedure outlined in equation (7). Other fracture and flow attributes are summarized in Table 5.Figure 12
Benchmark model, with (a) fracture geometry. Permeability ratios areRk=104. (b) Flow paths generated with EDFM and uniform pressure at the left-hand boundary is normalized with the right-hand boundary held at zero pressure. After [56].
(a)(b)Table 5
CAM model inputs.
Fracture/reservoir attributesSymbolValueNatural fracture width (m)W0.0002Natural fracture length (m)LGiven below∗Natural fracture height (m)H1Natural fracture angle to far-field flowα0°Porosityn0.1Effective far-field flow rate (m/s)Vx1.25×10−7Effective natural fracture strength (m4/s)υfScaled from equation (7)∗∗Matrix permeability (mD)k10Viscosity (cP)μ1∗Fracturelengths=15,26,30,33,35,51,64,69–71,76m. ∗∗Effectivenaturalfracturestrength=3.8,8.9,1.9,8.7,3.3,8.1,9.6,4.2,6.5,4.4,3.8,9.0,5.1×10−5m4/s.The results for the particle paths and pressure contours generated by the CAM fracture model are given in Figures13(a) and 13(b), respectively. The particle paths from CAM (Figure 13(a)) show acceptable, visual similarity to those generated with EDFM (Figure 12(b)). It should be emphasized that Shah et al. [56] evidently introduced new particle seeds for flow tracking behind certain fractures, which explains the increased density of streamlines in the right-hand half of the simulation area of Figure 12(b). Although we could introduce new particle seeds in CAM in any location, we have not done so in Figure 13(a). Any further differences between the flow paths of Figures 12(b) and 13(a) may be attributed to different boundary conditions at the top and bottom of the flow space. Additionally, some assumed variables such as fracture aperture may account for the local difference in the assigned permeability contrast (from equation (7)). Figure 13(b) shows the pressure contours (in Pa) generated from the CAM code.Figure 13
a) Particle paths for the model in Figure12 using CAM based on the inputs from Table 5. (b) Pressure contours (in Pa) for the model in (a) from CAM.
(a)(b)Also, trying to match permeability ratio scaling with numbers from the numerical, dual-porosity model (EDFM) may be misleading because the permeability ratio (Rk) stated as 104 in Shah et al. [56] may in fact be somewhat of a mixed input number. The permeability of the matrix is based on estimations conforming to Darcy’s law, but the permeability in the (nonporous) open fractures is based on a cubic equation (fracture width3) to obtain the right dimensions for the flux calculation in the open fracture. Therefore, our CAM model is fundamentally different from EDFM because we scale the permeability of the fracture based on an extension of the Darcy flow assumption, rather than using a cubic law (as for the open space fractures in EDFM).


### 5.2. CAM Runtimes

The steady-state pressure, velocity, and streamline contours can be computed instantaneously in CAM. However, the particle paths in CAM are generated from Eulerian particle tracking, which may be computationally intensive. The generalized Eulerian particle tracking algorithm is given by equation (8) [34]:
(8)zn+1t1≈znt0+Vznt0Δt,where zn is the initial position of a particle at time t0 with the velocity of Vzn. The particle paths are generated by first choosing an initial position z0 at time t0=0 and calculating the initial instantaneous velocity. By choosing an appropriate timestep, Δt, the position zjtj of the tracer at time tj is:
(9)zjtj≈zj−1tj−1+Vzj−1tj−1Δt.The runtimes (tcpu) for an Eulerian particle tracking scheme depend on the time step chosen for the simulation (Δt), the number of natural fractures present (Nf), the number of flow particles tracked (np), and the duration of simulated flow (tTotal) as summarized in equation (10).
(10)tcpu=fΔt,Nf,np,tTotal.The velocity (v) is assumed to be constant for Δt, which is a valid assumption for slow-moving fluids. The smaller time steps increase the accuracy of this discrete time approximation. When Δt is not small enough, the particles will overshoot the actual path to an adjoining path, closer or further from the original path [6]. The Δt must be selected by a trial and error approach, where a time step of unit time (e.g., 1 day) is initially chosen and then is reduced if particle paths are not smooth. The Δt needs to be reduced for stronger fractures to generate smooth particle paths. The next variable that affects the runtime of Eulerian particle tracking is the number of natural fractures (Nf) simulated. When the number of natural fractures increases, the velocity of the additional natural fracture needs to be superimposed, leading to a longer runtime. Next, the number of particles tracked also increases the runtime of the Eulerian particle tracking. The time-of-flight contours (TOFCs), which can be used to calculate the drained rock volume (DRV), are computed by connecting all the particle positions after a certain time period since the onset of flow. A densely seeded number of particles may be required, depending on the distance and strength of the individual natural fractures, to generate smooth TOFCs. Finally, the total run time depends on the duration of the flow simulation (tTotal), which determines the number of timesteps needed to complete the simulation.All the variables involved in Eulerian particle tracking need to be carefully selected on a case-to-case basis to optimize the speed and accuracy of the simulation. We investigated variousΔt for the model in Figure 13(a) to calculate the runtimes (tcpu), which is presented in Table 6 (MATLAB 2018b code on a Quad-core 3.4 GHz Intel i5-4670K).Table 6
Computational time (tcpu) for several test cases.
TimestepΔtdayNumber of fracturesNfNumber of particles trackednpDuration of simulationtTotalyearsTotal run timetcpus1136013042.60.5136013081.20.0513601307100.0213601301820.6


## 5.1. Model Description and Results

The model design and particle paths from the EDFM reference solutions of Shah et al. [56] were used as a starting point for comparison with CAM results (Figures 12(a) and 12(b)). The EDFM representation originated from a dual continuum model based on a square of unit dimensions transected by fractures with unstated apertures, and geometry as portrayed in Figure 12(a). The fracture-matrix permeability ratio Rk=kf/km which is constant for all fractures is 104. Figure 12(b) shows the flow paths generated from EDFM based on the boundary conditions in [56]. The boundary conditions in the benchmark of Figure 12(b) are a combination of two Neumann’s conditions (namely no-flow boundaries, such that the directional derivative vanishes) and two Dirichlet’s conditions (with directional derivatives normal to the left and right boundary being due to uniform pressures). Essentially, no-flow boundaries are imposed at the upper and lower boundaries of the 2D flow area studied, and constant, uniform pressures are maintained at the left and right boundaries.In our CAM model, a no-flow upper boundary and lower boundary can be simulated by a fracture channel, but this would require Schwarz-Christoffel’s and Schottky-Kleine’s prime function boundary mapping as was used in our models of bounded reservoirs [71]. However, CAM is most user-friendly when applied to unbound flow domains, which is why we assume a uniform flux entering the flow space of the unit square from the left boundary (akin to a uniform pressure) (Figure 12(b)). The fracture lengths and orientations are extracted from Figure 12(a). The width of the fractures is not given [56] and is assumed to be 0.2 mm (0.0002 m). Although the available field data for the width (or aperture) of a natural fracture are limited from subsurface observations, values of 0.01 to 10 mm are reported in the literature [64]. The permeability ratios are determined by scaling the fracture flow strengths following the procedure outlined in equation (7). Other fracture and flow attributes are summarized in Table 5.Figure 12
Benchmark model, with (a) fracture geometry. Permeability ratios areRk=104. (b) Flow paths generated with EDFM and uniform pressure at the left-hand boundary is normalized with the right-hand boundary held at zero pressure. After [56].
(a)(b)Table 5
CAM model inputs.
Fracture/reservoir attributesSymbolValueNatural fracture width (m)W0.0002Natural fracture length (m)LGiven below∗Natural fracture height (m)H1Natural fracture angle to far-field flowα0°Porosityn0.1Effective far-field flow rate (m/s)Vx1.25×10−7Effective natural fracture strength (m4/s)υfScaled from equation (7)∗∗Matrix permeability (mD)k10Viscosity (cP)μ1∗Fracturelengths=15,26,30,33,35,51,64,69–71,76m. ∗∗Effectivenaturalfracturestrength=3.8,8.9,1.9,8.7,3.3,8.1,9.6,4.2,6.5,4.4,3.8,9.0,5.1×10−5m4/s.The results for the particle paths and pressure contours generated by the CAM fracture model are given in Figures13(a) and 13(b), respectively. The particle paths from CAM (Figure 13(a)) show acceptable, visual similarity to those generated with EDFM (Figure 12(b)). It should be emphasized that Shah et al. [56] evidently introduced new particle seeds for flow tracking behind certain fractures, which explains the increased density of streamlines in the right-hand half of the simulation area of Figure 12(b). Although we could introduce new particle seeds in CAM in any location, we have not done so in Figure 13(a). Any further differences between the flow paths of Figures 12(b) and 13(a) may be attributed to different boundary conditions at the top and bottom of the flow space. Additionally, some assumed variables such as fracture aperture may account for the local difference in the assigned permeability contrast (from equation (7)). Figure 13(b) shows the pressure contours (in Pa) generated from the CAM code.Figure 13
a) Particle paths for the model in Figure12 using CAM based on the inputs from Table 5. (b) Pressure contours (in Pa) for the model in (a) from CAM.
(a)(b)Also, trying to match permeability ratio scaling with numbers from the numerical, dual-porosity model (EDFM) may be misleading because the permeability ratio (Rk) stated as 104 in Shah et al. [56] may in fact be somewhat of a mixed input number. The permeability of the matrix is based on estimations conforming to Darcy’s law, but the permeability in the (nonporous) open fractures is based on a cubic equation (fracture width3) to obtain the right dimensions for the flux calculation in the open fracture. Therefore, our CAM model is fundamentally different from EDFM because we scale the permeability of the fracture based on an extension of the Darcy flow assumption, rather than using a cubic law (as for the open space fractures in EDFM).


## 5.2. CAM Runtimes

The steady-state pressure, velocity, and streamline contours can be computed instantaneously in CAM. However, the particle paths in CAM are generated from Eulerian particle tracking, which may be computationally intensive. The generalized Eulerian particle tracking algorithm is given by equation (8) [34]:
(8)zn+1t1≈znt0+Vznt0Δt,where zn is the initial position of a particle at time t0 with the velocity of Vzn. The particle paths are generated by first choosing an initial position z0 at time t0=0 and calculating the initial instantaneous velocity. By choosing an appropriate timestep, Δt, the position zjtj of the tracer at time tj is:
(9)zjtj≈zj−1tj−1+Vzj−1tj−1Δt.The runtimes (tcpu) for an Eulerian particle tracking scheme depend on the time step chosen for the simulation (Δt), the number of natural fractures present (Nf), the number of flow particles tracked (np), and the duration of simulated flow (tTotal) as summarized in equation (10).
(10)tcpu=fΔt,Nf,np,tTotal.The velocity (v) is assumed to be constant for Δt, which is a valid assumption for slow-moving fluids. The smaller time steps increase the accuracy of this discrete time approximation. When Δt is not small enough, the particles will overshoot the actual path to an adjoining path, closer or further from the original path [6]. The Δt must be selected by a trial and error approach, where a time step of unit time (e.g., 1 day) is initially chosen and then is reduced if particle paths are not smooth. The Δt needs to be reduced for stronger fractures to generate smooth particle paths. The next variable that affects the runtime of Eulerian particle tracking is the number of natural fractures (Nf) simulated. When the number of natural fractures increases, the velocity of the additional natural fracture needs to be superimposed, leading to a longer runtime. Next, the number of particles tracked also increases the runtime of the Eulerian particle tracking. The time-of-flight contours (TOFCs), which can be used to calculate the drained rock volume (DRV), are computed by connecting all the particle positions after a certain time period since the onset of flow. A densely seeded number of particles may be required, depending on the distance and strength of the individual natural fractures, to generate smooth TOFCs. Finally, the total run time depends on the duration of the flow simulation (tTotal), which determines the number of timesteps needed to complete the simulation.All the variables involved in Eulerian particle tracking need to be carefully selected on a case-to-case basis to optimize the speed and accuracy of the simulation. We investigated variousΔt for the model in Figure 13(a) to calculate the runtimes (tcpu), which is presented in Table 6 (MATLAB 2018b code on a Quad-core 3.4 GHz Intel i5-4670K).Table 6
Computational time (tcpu) for several test cases.
TimestepΔtdayNumber of fracturesNfNumber of particles trackednpDuration of simulationtTotalyearsTotal run timetcpus1136013042.60.5136013081.20.0513601307100.0213601301820.6


## 6. Discussion

The current study provides an improved scaling rule (equation (7)) to model the natural fracture strength in CAM models when the permeability contrast is known. This study also includes a comparison of CAM model results with EDFM to complement our earlier validation [43]. Modern model efforts have exclusively employed finite-element methods due to certain perceived limitations of closed-form solutions. We claim that this narrow focus is unwarranted. The development of appropriately tailored closed-form solutions based on complex analysis methods (CAM) offers a number of unique strengths. These strengths include (1) highly compact formulation, (2) infinite resolution, and (3) ultrafast computation time (see Section 5.2). The infinite resolution is due to the lack of any gridding. Minor drift or dispersion may occur when the time-of-flight option is used, due to finite time-stepping increments. However, small timesteps will make any drift negligible. For the instantaneous pressure field and streamline solutions using the integral method, such effects do not occur.Some key differences in model design and associated input parameters of CAM and EDFM have come to the foreground in our study, which merits an in-depth discussion in Section6.1. A brief discussion about the fundamental assumption of irrotational flow, which should not be equated to a requirement of inviscid fluid flow, is given in Section 6.2. The potential merits of merging gridless solutions for flow and fracture propagation, when studying fractured porous media, are highlighted in Section 6.3.


### 6.1. Key Differences in Model Design and Parameters of CAM and EDFM

This article made a first attempt to compare CAM models for naturally fractured reservoirs with the results of EDFM. Several fundamental differences in model method design, which also reflect on input parameters, make a straightforward quantitative benchmark challenging. What we have offered in the present study is a qualitative comparison. The critical issues we run into when comparing CAM and EDFM are as follows.(1)
Natural fractures in CAM models are scaled in a fundamentally different way from those in EDFM(2)
EDFM is based on a dual continuum model assumption, which assumes that the fractures are open spaces with a pseudopermeability assigned with a value based on the fracture aperture using a cubic equation. Such open fractures are likely to have assigned permeability, being several orders of magnitude larger than the primary permeability(3)
CAM models scale the permeability of the fracture based on an extension of the Darcy flow assumption using a flux strength. The flux strength controls the amount of fluid that travels via the natural fractures, which can be varied also with the fracture aperture, but in CAM is not governed by a cubic law as in the EDFM model approach(4)
EDFM models commonly impose the Dirichlet boundary conditions with constant pressure, which in CAM is challenging to maintain, because the method works with a certain velocity and flux inputs with pressure being a consequence rather than a primary input parameter(5)
When natural fractures occur sufficiently remote from the Dirichlet boundaries, such as in Figures10 (EDFM) and Figure 13 (CAM), then the constant pressure and constant flux at the boundaries are nearly equivalent boundary conditions. Nonetheless, some differences exist, which only diminish when the flow distortions due to natural fracture systems are placed very remotely from the Dirichlet boundariesThe above observations explain why a quantitative benchmark of CAM models with natural fractures with independent discrete element-based methods (EDFM), involving differences in fundamental model assumptions (cubic law for fractures, finite boundaries, and others), will remain challenging. The present study is a first attempt to identify those challenges, such that future studies can refer to those differences and possibly come up with mitigating solutions.


### 6.2. Inviscid versus Viscous Flow

Fundamental arguments about why the application of potential theory should not be restricted to inviscid fluids have been highlighted by Weijermars [34]. The potential theory is also valid for viscous fluids (in addition to inviscid fluids) when the boundary-layer effects are minimal [34]. The analytical description of boundary-layer effects during the flow of fluids has been highlighted by Wang [72]. CAM models have been used to study the flow of terrestrial lava flows [34], with arguments given for relaxation of the inviscid constraint. However, many more arguments for potential theory not being limited strictly to inviscid flows were given by Joseph et al. [73]: “It is never necessary and typically not useful to put the viscosity of fluids in potential (irrotational) flow to zero.” This is cited to demonstrate the applicability of the potential theory to describe viscous flow, which is also supported by many detailed studies cited in [73].


### 6.3. Merging Gridless Solutions for Flow and Fracture Propagation

Recent efforts have shown that fracture propagation in elastic and poroelastic media can be modeled without resorting to finite difference or finite-element methods [74–76]. The development of the so-called time-stepped linear superposition method (TLSM) was motivated by the same gains underpinning CAM solutions for fluid flow: skipping of tedious gridding to achieve faster computation times, while preserving infinite resolution. TLSM allows the determination of nonplanar fracture propagation paths for multiple hydraulic fractures growing from the perforations during fracture treatment [74–76]. The proprietary, gridless, fracture propagation simulator (TLSM) can be coupled with the gridless simulator for fluid flow in fractured porous media (CAM). The fundamental theory and key algorithms for both types of simulators have been published in leading applied mathematics journals [17, 77–79].


## 6.1. Key Differences in Model Design and Parameters of CAM and EDFM

This article made a first attempt to compare CAM models for naturally fractured reservoirs with the results of EDFM. Several fundamental differences in model method design, which also reflect on input parameters, make a straightforward quantitative benchmark challenging. What we have offered in the present study is a qualitative comparison. The critical issues we run into when comparing CAM and EDFM are as follows.(1)
Natural fractures in CAM models are scaled in a fundamentally different way from those in EDFM(2)
EDFM is based on a dual continuum model assumption, which assumes that the fractures are open spaces with a pseudopermeability assigned with a value based on the fracture aperture using a cubic equation. Such open fractures are likely to have assigned permeability, being several orders of magnitude larger than the primary permeability(3)
CAM models scale the permeability of the fracture based on an extension of the Darcy flow assumption using a flux strength. The flux strength controls the amount of fluid that travels via the natural fractures, which can be varied also with the fracture aperture, but in CAM is not governed by a cubic law as in the EDFM model approach(4)
EDFM models commonly impose the Dirichlet boundary conditions with constant pressure, which in CAM is challenging to maintain, because the method works with a certain velocity and flux inputs with pressure being a consequence rather than a primary input parameter(5)
When natural fractures occur sufficiently remote from the Dirichlet boundaries, such as in Figures10 (EDFM) and Figure 13 (CAM), then the constant pressure and constant flux at the boundaries are nearly equivalent boundary conditions. Nonetheless, some differences exist, which only diminish when the flow distortions due to natural fracture systems are placed very remotely from the Dirichlet boundariesThe above observations explain why a quantitative benchmark of CAM models with natural fractures with independent discrete element-based methods (EDFM), involving differences in fundamental model assumptions (cubic law for fractures, finite boundaries, and others), will remain challenging. The present study is a first attempt to identify those challenges, such that future studies can refer to those differences and possibly come up with mitigating solutions.


## 6.2. Inviscid versus Viscous Flow

Fundamental arguments about why the application of potential theory should not be restricted to inviscid fluids have been highlighted by Weijermars [34]. The potential theory is also valid for viscous fluids (in addition to inviscid fluids) when the boundary-layer effects are minimal [34]. The analytical description of boundary-layer effects during the flow of fluids has been highlighted by Wang [72]. CAM models have been used to study the flow of terrestrial lava flows [34], with arguments given for relaxation of the inviscid constraint. However, many more arguments for potential theory not being limited strictly to inviscid flows were given by Joseph et al. [73]: “It is never necessary and typically not useful to put the viscosity of fluids in potential (irrotational) flow to zero.” This is cited to demonstrate the applicability of the potential theory to describe viscous flow, which is also supported by many detailed studies cited in [73].


## 6.3. Merging Gridless Solutions for Flow and Fracture Propagation

Recent efforts have shown that fracture propagation in elastic and poroelastic media can be modeled without resorting to finite difference or finite-element methods [74–76]. The development of the so-called time-stepped linear superposition method (TLSM) was motivated by the same gains underpinning CAM solutions for fluid flow: skipping of tedious gridding to achieve faster computation times, while preserving infinite resolution. TLSM allows the determination of nonplanar fracture propagation paths for multiple hydraulic fractures growing from the perforations during fracture treatment [74–76]. The proprietary, gridless, fracture propagation simulator (TLSM) can be coupled with the gridless simulator for fluid flow in fractured porous media (CAM). The fundamental theory and key algorithms for both types of simulators have been published in leading applied mathematics journals [17, 77–79].


## 7. Conclusions

In this study, we presented a correction factor for scaling the permeability contrast of discrete natural fracture elements in CAM models, when the permeability contrast for a naturally fractured reservoir is known. Previously, the ratio of the maximum fluid velocity in natural fractures relative to the unperturbed far-field velocity was understood to be a proxy for permeability contrast in CAM models. A comparison with a numerical model showed that the ratio of fluid velocity underestimated the actual permeability contrast by a factor 2.5 and needs a correction factor. The particle paths generated from CAM were compared with similar results from EDFM. In addition to scaling issues, this study investigated the effect of fracture aperture and permeability changes on pressure contour plots. The following conclusions can be drawn from our study:(a)
CAM can be used to model the flow of fluids in fractured porous media. The natural fractures locally increase fluid velocity or decrease it according to the perviousness of the natural fractures. For a conductive fracture, the flow outside of the matrix remains largely unaffected (Figure2(b)), although particle tracking shows that considerable distortion occurs locally (Figures 6(b) and 9). The impervious fracture (Figure 10) also distorts the pressure contours and changes the path of fluid flow in the vicinity of the fracture. The TOFCs in the wake of the cemented fractures are pulled back due to the fluid particles being slowed down(b)
CAM models allow high-resolution visualization of particle paths, pressure, and velocity fields without complex gridding and meshing. The presence of natural fractures in a reservoir promotes preferential flow channeling (Figure13). This can help in faster transport of fluid from the matrix to the wellbore. However, it can also increase fracture and well interference due to pressure communication between the closely spaced wells and hydraulic fractures(c)
The permeability contrast calculated from the ratio of matrix and fracture permeability needs a correction factor of 2.5 to calculate the strength variable for CAM(d)
The particle paths generated from CAM were compared to the results from EDFM. Despite of differences in boundary conditions and other assumptions, the results from both models showed a visual match. Quantitative benchmarks of the CAM model of naturally fractured reservoirs with EDFM or other discrete volume methods will remain challenging because of the fundamental differences in the design assumptions (explained in Section6.1)



---

*Source: 8838540-2020-10-13.xml*