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

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.


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.

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 lowpermeability and low-porosity fractured media. For a lowporosity/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 2 Geofluids 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.

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][36][37][38]. Models based on CAM have been previously used to model fluid flow in hydraulically fractured reservoirs [39][40][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.

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 3 Geofluids 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 m 2 /s, the strength of a hydraulic fracture has units of m 3 /s, and the strength of an areal doublet/natural fracture element has units of m 4 /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].

CAM Strength Scaling of Natural Fractures. From
Darcy's law, time-dependent flow rate, V f , across a natural fracture of length L, due to the time-dependent pressure gradient ΔP f ðtÞ is defined as follows [17]: where k f /μ is the ratio of fracture permeability to fluid viscosity. Similarly, the time-dependent flow rate, V m , across a section of the matrix with the same length due to the time-dependent pressure gradient ΔP m ðtÞ is defined as follows [17]: Assuming that the pressure gradient across the natural fracture and adjoining matrix is equal (i.e., ΔP m ðtÞ = ΔP f ðtÞ), from equations (1) and (2), the permeability ratio k f : k m = R k can be calculated as follows [17]: 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 R k of the matrix and the natural fracture can be calculated from the ratio of the respective fluid velocities.
In our previous study, the permeability ratio R k was linked to the primary input parameters of CAM [6], which are modified here to account for the superposition of V f being the result of V m and the superposed flux. Our prior study [6] where V m is the far-field velocity in the matrix, and H f , W f , and L f are the height, width, and length of the natural fracture element, respectively. However, the velocity in the fracture, V f , is due to a preexisting V m plus an additional velocity component V flux , superposed due to the fracture flux: Substituting now the modified fracture strength (4) and using R k yields the following: The strength of the natural fracture and the far-field velocity in the matrix are denoted by υ f ′ and V m , 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: 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 (R k ) calculated from equation (3) underestimates the permeability contrast 4 Geofluids 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 (k f /k m ) 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.

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 of 1:117 × 10 −7 m/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 −6 m/s. In the remainder of this study, the effective strength and the effective velocity are reported, accounting for the porosity of the reservoir. Next, consider a reservoir model ( Figure 2) 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.
Figures 2(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 (R k ) 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 Figure 2(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 ( Figure 2(d)) is the primary reason for preferential flow paths and flow channeling due to fracture networks [49][50][51][52]. Multiple tracer transport studies on core samples of different length scales have shown that flow inside a fractured reser-voir is highly heterogeneous [52][53][54][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]. Figure 3 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.  Table 3 (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

Geofluids
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.
EDFM is a special form of a discrete fracture model (DFM) model (see Section 2.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]. Figure 4 considers a single-phase flow in a 9 × 9 m 2 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 m 2 . The flux for each (matrix-matrix and fracture-fracture) grid interface is defined by using the twopoint-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,   [56][57][58][59].
The flow in Figure 4 is driven by the Dirichlet boundary conditions of 10 7 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 (R k ) 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 ( Figure 4) 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  (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 R k for the model generated from EDFM (Figure 4(a)) is 50, whereas the R k for CAM is 20 ( Figure 5(d)). Hence, an empirical correction factor ξ is introduced to scale the R k 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: where R k 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).

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.      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.
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   (7), where ξ = 2:5. 8 Geofluids aperture sizes. Section 4.2 presents the pressure contours for natural fractures with reduced permeability relative to the ambient matrix rock in the reservoir.

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][64][65]. In this section, we investigate the effect of fracture width and orientation on pressure contour patterns using a single natural fracture. Figure 7 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 Table 4. The pressure contour patterns for both Cases A and B are distorted near the fracture tips (Figures 7(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). The particle paths represent the progressive movement of fluid over 30 years. The particle paths (blue) and the timeof-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 R k constant but adjust the flux strength υ f in proportion to the fracture width W f , 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 R k . A larger R k will promote fracture channeling and lead to a shorter TOFC. When R k < 1, the fracture becomes progressively impervious, leading to a longer TOFC for fluid traveling via the fracture. Also, when R k = 1, the fracture may physically exist, but its presence will not affect the flow paths.

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  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 −6 m 4 /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). 9 Geofluids natural fracture in Case A (Figure 7(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 −7 m 4 /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.

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  10 Geofluids 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   11 Geofluids 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).

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 R k = k f / k m which is constant for all fractures is 10 4 . Figure 12(b) shows  12 Geofluids 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, noflow 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 userfriendly 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.
The results for the particle paths and pressure contours generated by the CAM fracture model are given in Figures 13(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 fur-ther 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.
Also, trying to match permeability ratio scaling with numbers from the numerical, dual-porosity model (EDFM) may be misleading because the permeability ratio (R k ) stated as 10 4 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 width 3 ) 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 (a) (b) Figure 12: Benchmark model, with (a) fracture geometry. Permeability ratios are R k = 10 4 . (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].

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]: where z n is the initial position of a particle at time t 0 with the velocity of Vðz n Þ . The particle paths are generated by first choosing an initial position z 0 at time t 0 = 0 and calculating the initial instantaneous velocity. By choosing an appropriate timestep, Δt, the position z j ðt j Þ of the tracer at time t j is: The runtimes (t cpu ) for an Eulerian particle tracking scheme depend on the time step chosen for the simulation (Δt), the number of natural fractures present (N f ), the number of flow particles tracked (n p ), and the duration of simulated flow (t Total ) as summarized in equation (10).
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 (N f ) 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 (t Total ), 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 (t cpu ), which is presented in Table 6 (MATLAB 2018b code on a Quad-core 3.4 GHz Intel i5-4670K).

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  Figure 13: a) Particle paths for the model in Figure 12 using CAM based on the inputs from 14 Geofluids (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 Section 6.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 The 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.

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 boundarylayer 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].

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][75][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: 15 Geofluids 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][75][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][78][79].

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 (Figure 2(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 ( Figure 13). 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 dif-ferences 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 Section 6.1) θ 1 = π/2 ð Þ, θ 2 = −θ 1 , γ 2 = γ 1 − π/2 ð Þ,

B. The Effect of Increased Fracture Strength on Branch Cuts
One crucial aspect of CAM models is the occurrence of mathematical branch cuts when multivalued solutions appear along certain integral lines [7,44]. Prior studies have discussed possible solutions to side-step such branch-cut effects [7,44]. When the fracture strength superposed on a far-field flow is increased, the appearance of pressure jumps across branch cuts signals that the flow becomes physically unrealistic. For instance, Figure 8 shows the pressure contours for fractures, where the effective fracture strength is doubled to 4:76 × 10 −6 m 4 /s (as compared to 2:38 × 10 −6 m 4 /s in Figure 7(a), main text). The branch-cut effect is nearly negligible for the case where the fracture is parallel to the direction of fluid flow (Figure 8(a)). However, the branch cuts cause a discontinuity in pressure contours even when the fractures are slanted (Figures 8(b) and 8(c)). CAM is an analytical model where the fracture strength may be set to physically unrealistic values. Figure 11 shows the pressure contours for a blocking fracture with a reduced effective fracture strength (−1:92 × 10 −6 m 4 /s). The pressure contours of Figure 11 show the same pattern as for the blocking fracture in Figure 9 (main text), but the effect of the branch cut is more apparent when the strength decreases, as shown in Figure 11.

Data Availability
The underlying data is included in the paper.

Conflicts of Interest
The authors declare that they have no conflict of interest.