Flow in Fractured Porous Media Modeled in Closed-Form: Augmentation of Prior Solution and Side-Stepping Inconvenient Branch Cut Locations

Carefully chosen complex variable formulations can solve flow in fractured porous media. Such a calculus approach is attractive, because the gridless method allows for fast, high-resolution model results. Previously developed complex potentials to describe flow in porous media with discrete heterogeneities such as natural fractures can be modified to expand the accuracy of the solution range. The prior solution became increasingly inaccurate for flows with fractures oriented at larger angles with respect to the far-field flow. The modified solution, presented here, based on complex analysis methods (CAM), removes the limitation of the earlier solution. Benefits of the CAM model are (1) infinite resolution, and (2) speed of use, as no gridding is required. Being gridless and meshless, the CAM model is computationally faster than integration methods based on solutions across discrete volumes. However, branch cut effects may occur in impractical locations due to mathematical singularities. This paper demonstrates how the augmented formulation corrects physically unfeasible refraction of streamlines across high-permeability bands (natural fractures) oriented at high angles with respect to a far-field flow. The current solution is an important repair. An application shows how a drained rock volume in hydraulically fractured hydrocarbon wells will be affected by the presence of natural fractures.


Introduction
Fast and accurate solutions for flow in fractured porous media are sought after in the petroleum industry, which increasingly extracts oil and gas from hydraulically and naturally fractured shale formations. Practical applications for flow models of fractured porous media are ubiquitous and not limited to the petroleum industry. For example, medical applications involving flow across fractured bone structures, and magmatic melts moving from the Earth's mantle to the crust via pore space and fractures also are governed by equations for flow in fractured porous media. Of particular interest to the current study is the modeling of the flow of hydrocarbon fluids in fractured reservoirs, for which our models can provide detailed solutions. However, the formulation provided here may also be applied to subsurface flow studies for geothermal energy extraction from naturally fractured rocks.
Modeling flow in fractured porous media with discrete fractures using explicit numerical methods is computationally intensive and requires excessive meshing and gridding adjustments near fracture intersections [1][2][3][4]]. An alternative model approach, advocated here, uses a gridless and meshless analytical modeling technique based on complex analysis methods (CAM), and has been stepwise developed in two key, earlier papers. The first paper [5] presented a solution using areal doublets,

Complex Potential Solutions
In this section, we revisit some fundamental closed-form solutions required to ultimately adjust complex potentials used to describe flow through discrete fracture systems. Redirection of streamlines can be effectuated when an existing flow, described by a set of single-valued complex potentials is superposed by a complex potential that represents the desired disturbance. For example, Figure 1a shows an undisturbed flow between a direct line drive with streamlines connecting a series of injectors and producer wells. The basic flow pattern will be altered when a rectangular domain of high permeability is placed in the flow path between the injector and producers; the streamlines converge (with a faster time-of-flight) toward the producer wells ( Figure 1b). The opposite occurs (i.e., a slower time-of-flight) when the rectangular domain instead has a lower permeability than the surrounding matrix, which will result in streamline divergence (Figure 1c). Both the flow convergence and divergence (in Figure 1b,c) were modeled by superposing a series of line doublets [7]. Section 2 moves through three systematic steps that highlight how the line doublets and related line dipoles are transformed into areal doublets and line dipoles that are suitable for modeling flow through discrete natural fractures in a reservoir. First, Section 2.1 explains the difference between the line doublet and a line dipole, which are two end members of a continuous series of line doublet/line dipole intermediates. Next, Section 2.2 reviews the prior complex potential used for closed-form flow models of fractured porous media and gives a comprehensive complex potential which can describe both the areal doublet and areal dipole. Section 2.3 highlights the limitations to the practical use of the derived solution due to branch cut effects. The results of Section 2 are necessary for an amended solution given in Section 3 that avoids the branch cut effects. (a-c) Waterfront advance between a series of injector outlets and producer wells in streamline visualizations (left column; yellow flow paths inside grey injection water body; blue curves beyond injection water body). Time-of-flight contours (TOFCs, red curves). Velocity field magnitude (right), scaled by a Bernoulli pressure proxy using the natural logarithm and inverse tangent, outlines the regions (blue) where the fastest flow occurs. Reproduced with permission from [7].

Line Doublets and Line Dipoles
The distinction between doublets and dipoles was originally proposed by Strack [8], and was required for constructing line arrays called, respectively, line doublets and line dipoles.   The complex potential for a line doublet and a line dipole are invariant [6], which means these potentials represent two physically distinct flows. The comprehensive solution presented below shows that the invariant nature of the respective complex potentials does not preclude the formulation of a concise complex potential that covers a broader spectrum of physical flows, including the line doublet and line dipole as end-members. The comprehensive solution is given by the following complex potential (Ω) [5]:  (1) is to introduce z a = z c − 0.5We iβ and z b = z c + 0.5We iβ [9], where W, z c and β are the width, center coordinates and tilt angle of the line doublet/dipole element. The need to have two complex coordinates (z a and z b ) is removed, and instead the mid-point (z c ) between the two ends of the region occupied by the flow element is given by specifying a width W. The polarity of the singularity doublet/dipole can be controlled by θ. For θ = 90 • , we obtain the line doublet solution given in Weijermars and Van Harmelen [9]. However, when β = 0 • , the angle θ of the doublet singularities in the line element does not need to be 90 • . For θ = 0 • we have the line dipole solution. Intermediate flow types occur for 90 • < θ < 0 • (see below). There are two separate methods to obtain the streamlines for flows described by the complex potential of Equation (1). The first method separates the real and imaginary parts of the complex potential to obtain, respectively, the potential function and the stream function. The solutions for the stream function and the potential function can then be obtained (for any steady-state flow moment at time t) as an integral solution with standard software packages (like MATLAB, Python, Mathematica or others), as follows:  Table 1, and various angles for doublet polarity given by θ.    The second method to plot the streamlines for the flow of Equation (1) is to trace individual particles with an Eulerian particle tracking method. The particle tracking for any particular flow is made possible by first deriving the velocity field from the complex potential specific to the flow of concern, using: Application of the derivative of Equation (3a) to the complex potential of Equation (1), using z = x + iy = z (cos θ + i sin θ) = r(cos θ + i sin θ) = re iθ , gives: Equation (3b) is the generalized line doublet/dipole equation, which can be used to generate the velocity field for any line doublet/dipole analytical element with an angle of β to the real axis measured in an anti-clockwise direction. Additionally, θ can be used to change the angle of the polarity of the singularities, which collectively makes the line element. As before, for a special case when θ = 90 • , e iθ = i, this equation becomes identical to Equation (C9) for the line doublet in Weijermars and van Harmelen [9] and Equation (B2) in van Harmelen and Weijermars [5]. Note that the complex conjugate z = re −iθ , and Euler's formula e iθ = cos θ + i sin θ gives for θ = 0 • e iθ = 1, and for θ= 90 • it follows that e iθ = i. Equation (3b) gives (v x , v y ) for fluid particles in each location of the flow field. The particle tracking method uses a time-step, ∆t, and initial particle position z 0 at time t 0 = 0 to track the streamlines followed by the particles. The position z 1 (t 1 ) of the tracer at time t 1 = t 0 + ∆t is: The position z j (t j ) of the tracer particles at any other future time t j is: Figure 4a-d show the streamlines for the complex potential of Equation (1), based on the particle tracking method of Equations (3b) and (4a,b) using selected angles θ and tilt angle of the line element β = 0 • . Streamline solutions of Equations (3b) using the Eulerian particle tracking method of Equations (4a,b) offer the advantage that so-called time-of-flight contours (TOFCs) can be constructed [7]. Additionally, the flow is integrated over time and the particle paths represent the true flow path of time-dependent flows. The TOFCs show the advancement of an infinite number of tracers for periodic (constant) time-steps.
Comparison of Figures 3a-d and 4a-d reveals that the two methods provide the same streamline solutions. The first approach (stream function approach; Figure 3a-d) gives additional information on the fluid fluxes (stream function values relative to the origin). The alternative particle tracking method (Figure 4a-d) gives additional information on the time-of-flight, which is useful for studies of hydrocarbon drainage and reservoir depletion [7], studies of geothermal reservoirs that require monitoring of the advancement of injected fluids [9], and for studies of subsurface flow including the tracking of pollution plumes.  A prior study has integrated an infinite number of line doublets to formulate a complex potential for an areal doublet ( Figure 5). The flow field of an areal doublet can be interpreted as an infinite number of doublet singularities with poles aligned with the direction of flow as schematically shown in Figure 5a. The complex potential Ω 1 (z, t) for the areal doublet, scaled by a flow strength υ(t) [m 4 s −1 ] is [5]:

Areal Doublets and Areal Dipoles
with dimensions of the areal doublet given by reservoir depth, H, doublet length, L, and aperture (or width), W, with porosity, n. The dummy terms A − D in the real and imaginary parts [ , for an areal doublet with corner points located at z a1 , z a2 , z b1 , and z b2 , and tilt angle γ with the real axis ( Figure 6), are as follows: Fluids 2020, 5, 51 With corner point co-ordinates given by: A generic, comprehensive complex potential can be formulated for both areal doublets/areal dipoles (Figure 5a,b), and areal flow elements with singularity doublets/dipoles oriented oblique to the boundaries of the areal flow element. In analogy to Equation (1), the prior solution of an areal doublet [5,6] given here in Equation (5), can be modified into a general complex potential for areal doublets/dipoles (Figure 5a), as follows: With the corner-point coordinates again given by Equation (7). The solution of Equation (8) differs from the prior solution of Equation (5) in that the doublet/dipole polarity is separately parameterized, by not necessarily letting θ = β. The velocity potential corresponding to Equation (8) is given by: This is the general equation that allows for setting the angle of all the analytical elements (point dipole/doublet (θ) to line dipole/doublet (β) to areal dipole/doublet (γ)). In Equations (6) and (7), we have the imaginary i in the equation and we have e iβ in the denominator, which will cancel each other out for all the cases, and is kept only to maintain analogy with the previous Equations (1) and (5).
The streamlines can again be obtained from Equation (8) (9), for a range of polarity angles: with a tilt angle for the line doublet element β = 90 • . Stream function values (ψ) given by rainbow colors (m 2 ·s −1 ) relative to line y = 0. Branch cuts appear at the fracture tips. Solved using the imaginary part of Equation (8).

Branch Cut Effects
Equation (8) correctly solves the streamlines for a generic areal doublet/dipole element and can be used to represent the streamline deflections occurring near high permeability fractures in fractured porous media. The solutions of Figure 8 show the progressive rotation of a fracture element and the stationary far-field flow (superposed on the doublet-dipole flow) by selecting θ such that θ + γ = 180 • , for a variety of angles γ. The composite complex potential, Ω Comp (z, t), including the complex potential for the uniform far-field flow, Ω 4 (z, t) = e iθ u ∞ (t)z (with u ∞ (t) already being the effective velocity, accounting for porosity) is: The streamline solutions given in Figures 7 and 8 are perturbed by the branch cuts emanating from the singularities, which cause "apparent" jumps in the stream function values from positive to negative, and vice versa. For the areal doublet (Figure 5a), the branch cut seems not to occur, but these are present as a row of singularities in the high and low streamline values. The corresponding potential function solutions also exhibit branch cuts jumps, which we discussed in considerable detail in a prior study [6].
The issue we now face is to move the branch cuts for the general solution to a place where they would not render the use of the Equations (8) and (9) impractical. For the pure areal doublet (Figure 5a), we have earlier succeeded to find a practical branch cut placement solution [6]. However, as pointed out by Holzbecher [10], available procedures for removing branch cut effects from analytical solutions remain limited.
In spite of numerous efforts in our present research program to move the branch cuts from the general solutions of Equations (8) and (9) to a practical location, we did not succeed, which meant that we could not apply Equations (8) and (9) in practical situations that required continuous Eulerian particle tracking. The next section discusses an alternative, augmented solution for areal doublet representation of flow deflections by discrete natural fractures of arbitrary orientation, which solves the prior limitation.

Augmented Solution
A new flow element used to model flow in fractured porous media was generated by the integration of line doublets [5]. The derived flow element was initially termed a natural fracture element [5], but later renamed as an areal doublet [6]. The reason for the re-examination of the earlier approach is that the existing formulation is accurate for fractures aligned with a far-field flow, and works well for fractures oriented moderately obliquely with respect to the far-field flow, but becomes increasingly inaccurate when the fractures are perpendicular to or at a large angle to the far-field flow. Figure 9 illustrates the problem arising when the simple area doublet formulation of Equation (5) is used to progressively increase the angle of the natural fracture (represented by a narrow areal doublet) relative to a steady far-field flow. A long and narrow fracture aligned with the far-field flow direction can be modeled accurately with an areal doublet ( Figure 9a). However, for the larger incidence angles, the flow solution is no longer accurate. Figure 9f shows the solution for a fracture at a large angle to the far-field. The areal doublet would require reorientation of the polarities in the fracture to align with the far-field flow. If we try to implement such polarity changes using Equation (9), the branch cut effects do not allow continuous particle tracking across the fracture region.

Shortcomings of the Prior Solution
To develop an augmented solution, while avoiding the branch cut effects, we first discuss the two end-members, areal doublet and areal dipoles, in some detail (Sections 3.2 and 3.3) paying specific attention to the orientation of corner point co-ordinates given in Equation (7). Ultimately, the branch cut issues that rendered Equation (9) impractical for use in real world problems were sidestepped by superposing the two end-members of the areal doublet/dipole family of solutions.

Areal Doublet Solution
The stream function solutions of the areal doublet of Equation (5) with the integral method are given in Figure 10a,b for the two extreme orientations of the natural fracture prototype (vertical and horizontal, as in Figure 9a,f), with a superposed far-field flow. Figure 10c,d show the areal doublet flow without the superposed far-field flow. The solutions given in Figure 10 are possible with Equation (5) while avoiding undue branch cuts by changing the corner point co-ordinates given in Equation (7), according to the respective L/W pairs of the fracture slots (as indicated in the sketches next to the respective solutions). When applied to Eulerian particle tracking, Equation (5) gives continuous solutions, as illustrated in the example of Figure 9. However, for the larger inclinations of the natural fracture element with respect to the far-field flow direction (e.g., Figure 9e,f) the solutions for the flow paths become unrealistic. The initial repair of the natural fracture algorithm occurred via the steps in Section 2.2. However, the proposed algorithms of Equation (8) resulted in multiple branch cuts, which make the continuous Eulerian tracking of particle movement impossible to use. We investigated numerous avenues to move branch cuts to locations where they would not cause multivalued functions along the particle paths, using our experience of a prior study [6]. All efforts failed, except for one approach via the steps explained in Sections 3.3 and 3.4.

Areal Dipole Solution
The alternative model for the long fracture aligned with the far-field flow is given by the complex potential for the areal dipole [6]: Equations (5) and (11) differ in that their real and imaginary axes are switched, which is achieved by multiplying the complex potential of Equation (5) with the imaginary number i. The dummy terms A − D remain unchanged and are given by Equations (6a)-(6d). Figure 11a-d are solutions for the areal dipole for the same orientations as in Figure 10a-d, with and without a superposed far-field flow. The solutions given in Figure 11 are possible with Equation (11) while avoiding undue branch cuts by changing the corner point co-ordinates given in Equation (7), according to the respective L/W pairs of the fracture slots (as indicated in the sketches next to the respective solutions). Although the solutions of Figures 10 and 11 appear the same, the important difference lies in the coordinate transformation, due to which the orientations of the examples shown in each set of Figures are swapped 90 • . This is a simple but essential step to correctly model the flow across natural fracture elements at a high angle to the far-field flow (see Section 3.4). The branch-cut effects that impeded continuous particle tracking across the natural fracture element with the solution of Equations (8) and (9)

Superpositions
Changing the corner point co-ordinates according to the respective L/W pairs of the fracture slots as applied in Figures 10 and 11 will avoid the occurrence of branch cuts in impractical places, but is not very user-friendly when superposing multiple fracture elements in practical applications. Thus, the velocity fields for two areal doublets are superposed with the velocity field of the far field flow, such that the two end angles (0 • and 90 • ) result in the correct particle-path solutions based on the conceptual model of Figure 12.
(−1) n+1 sin γ n ·Ω n (z) (12) Ω f (z) is complex potential for the far-field flow given by: Ω n (z) is the complex potential for the areal doublet/dipole element given by: The vertices z a1n , z a2n , z b1n , and z b2n are the vertices for the areal doublet/dipole element given by: The transformations and inputs needed to calculate the vertices given in Equation (15) are given below: Figure 13 shows the particle path across a single, highly conductive fracture at different orientations with respect to a far-field flow incoming from below generated by using the velocity potential derived from Equation (12). The length (L 1 ) and the width (W 1 ) were assumed to be 10 m and 1 m respectively. The results in Figure 13 show more plausible flow paths compared to the results in Figure 9 for the cases where the fracture is tilted at an angle of more than 45 • with respect to the direction of the far-field flow.

Application
In this section, the augmented solution for flow across discrete fractures will be applied to a practical situation. Natural fractures may have a significant impact on the performance and flow of hydraulically fractured wells tapping into hydrocarbon reservoirs [11]. The interaction of hydraulic and natural fractures makes forecasting the performance of wells in such reservoirs a major challenge [12,13]. Alternatively, conductive natural fractures may promote pressure communication between the adjoining hydraulic fractures and wells, which increase the well/fracture interference [14]. The nonconductive natural fractures may affect the sweep pattern due to waterflood between the injector and the producer wells [15]. Hence, the estimation of a reservoir's flow performance and overall production depends on accurate description of the flow diversion due to the presence of natural fractures and natural fracture networks. Section 4.1 provides a brief background on different modeling methods for naturally fractured reservoirs.
Natural fractures distort the drained rock volume (DRV) around hydraulic fractures and may result in pressure communication with fractures of the adjoining stages and wells. The DRV is defined as the reservoir volume that contributes to oil and gas production through the movement of fluid particles to the wellbore after a particular specified production period. Although there is movement of fluid particles throughout the reservoir, only a limited fraction of those particles reaches the wellbore in a reasonable time period, which is shown by the DRV. Unusual refractions of streamlines near the well system occurred in our prior model code (compare Figures 9f and 13f) when natural fractures were oriented at a high angle to the far-field flow direction. We revisit the examples presented in an earlier study (in Section 4.2), where the effect of natural fractures on the DRV was investigated by using the algorithm presented in Equation (5) [16]. In Section 4.3, the DRV is regenerated using the modified solution of Equation (12). Finally, in Section 4.4, we compare results obtained earlier [16] with the results from the modified algorithm presented in this paper. The DRV for the reservoirs where natural fractures are oriented in the direction of flow remains unchanged. However, the DRV for the reservoirs where the natural fractures are oriented at a high angle to the far-field flow direction show a drastic change. A brief validation of the augmented solution is given in Section 4.5.

Flow in Fractured Reservoirs
Oil and gas reservoirs may contain vast networks of natural fractures, which interact with the hydraulic fractures, and therefore the most accurate flow simulation and production forecasting models must strive to account for the impact of natural fractures, based on their orientation, distribution, connectivity, strength and interaction with the hydraulic fractures [17,18]. Previously, naturally fractured reservoirs were mainly modeled by using dual porosity and modifications of the dual porosity model [19,20]. These models consider natural fractures as a secondary continuum with up-scaled flow properties, but have limited accuracy, and cannot be used to model discrete fracture networks [21]. Another widely used technique uses discrete fracture network (DFN) models [22] and embedded discrete fracture network models (EDFM) [23]. Such discrete models have some benefits over dual porosity models, such as flexibility and give more accurate representations of the fracture network, but they are computationally very demanding and expensive for complex reservoirs with numerous natural fractures. Both the continuum model (dual porosity) and the discrete model (DFN, EDFM) have been thoroughly reviewed in our prior study [16].
Hence, there is a need for new fast and gridless methods (such as the one presented in this study), which allow for efficient modeling of fluid flow in naturally fractured reservoirs. The present study builds forward on the earlier studies where closed-form analytical solutions based on complex analysis methods (CAM) were applied to single-phase flow reservoirs, assuming either hydraulic fractures only [24] or hydraulic fractures intersecting with natural fractures [6]. The natural fracture algorithms have been used in fundamental studies of the equivalent tensor concept [25]. An advantage of CAM models over numerical models is the ability to visualize the drained rock volume (DRV) with infinite resolution.
Unlike discretized numerical simulators, the CAM simulator offers enhanced visual resolution of the DRV around the hydraulic and natural fractures. Flow visualization at high resolution can be especially important to optimize fracture and well spacing in already crowded shale fields with several horizontal wells and hydraulic fractures. The prior studies [24,25] used an algorithm developed earlier [5], which works well when the natural fractures are oriented in the direction of the fluid flow. However, the algorithm needs to be modified when the natural fractures occur at larger angles to the direction of the fluid flow, as illustrated in the following section.

Field Application: Flow Near Hydraulic Fractures
In this section, we revisit some results of our previous study [16], which presented an extensive collection of results where natural fractures of varying orientations were introduced to a base case (Figure 14), and were simulated using Equation (5).  Figure 14a shows the streamlines (blue color) for a full-scale model of a hydraulically fractured well. The particle paths are generated only for the middle five hydraulic fractures for computational efficiency. However, the flow toward all the hydraulic fractures in the complex plane (119 fractures) is accounted for in the model simulation and were simulated for a period of 30 years. The identical hydraulic fractures, each of which has a fracture half-length of 150 ft (45.7 m) and height of 60 ft (18.3 m), are assumed to have a constant flux of 13.6 ft 3 /day (0.39 m 3 /day) per fracture, which is equivalent to the hydraulic fracture (line sink) strength of 8.47 ft 2 /day (0.79 m 2 /day). The rainbow colors in Figure 14b show the time of flight contours (TOFC), where each different color represents the DRV expansion for three-year increments from five hydraulic fractures (one stage) in the center of the well.
When natural fractures are present in the reservoir, the DRV shape around the hydraulic fractures will be progressively distorted when natural fractures (oriented oblique-at 45 • -to the horizontal well) occur in the proximity of the hydraulic fractures. The permeability contrast between the matrix and the fractures is increasing from the left to the right ( Table 2). Figure 15c,d shows the DRV distortion due to highly conductive natural fractures oriented at a high angle with respect to the far-field flow (10 • ) to the horizontal well). The natural fracture permeability contrast with the matrix (and corresponding strengths) and other properties for each of the cases are presented in Table 2. The DRV, for each case (Figure 15a-d), is distorted. However, where natural fractures have steep angles with respect to the far-field flow, the resulting sideward deflection of the particle paths appeared physically unrealistic (Figure 15c,d). As shown earlier (in Figure 9), the application of Equation (5) leads to such unrealistic particle paths for steep natural fractures. For higher angles to the far-field flow direction, the results of Figure 15c,d present an improved solution (Section 4.3).  Figures 14 and 15. Modified from [8].

Augmented Solution for Flow Near Hydraulic Fractures with Natural Fractures
Next, we implement the augmented solution of Equation (12) and rerun the results, aimed at removing any unrealistic particle paths seen in Figure 15. Figure 16a-d shows the DRV generated from the augmented solutions corresponding to the cases shown in Figure 15a-d. Comparison of the DRVs for the first two cases (Figure 15a,b and Figure 16a,b) show minimal effect due to the change in algorithms because the natural fractures are oriented at a moderate angle of 45 • to the direction of the particle paths. However, the DRV distortion intensifies for highly conductive natural fractures oriented at a large angle to the far-field flow direction (Figure 15c,d), as can be calculated from the augmented equation (Equation (12)). Table 2. Input summary for reservoir and natural fractures simulated in Figures 14 and 15. Modified from [8].

Comparison of Results Augmented and Earlier Solution
The DRV shapes of Figure 15c,d and Figure 16c,d differ significantly. The extreme deflection of particle paths seen in Figure 15c,d does not occur in Figure 16c,d, due to the application of the augmented code based on Equation (12). The extent of the fluid plume drained from the reservoir reaches farther into the reservoir (up to 1000 ft or 305 m, with natural fractures) for all cases (even after the application of the augmented solution), as compared to only 800 ft (244 m) in a reservoir without natural fractures (Figure 14b) or with few natural fractures oriented at moderate angles to the far-field flow (Figure 15a,b and Figure 16a,b).

Verification
Finally, the augmented solution is verified here in a straightforward experiment, which compares the prior solution of Equation (5) with the augmented solution of Equation (12). Consider a square domain with a higher permeability than the ambient matrix. The flow through the domain can be simulated by superposing an areal doublet of a certain strength on the far-field flow (Table 3), similar to what is shown in Figure 1b. The left column in Figure 17 gives the prior solution for three angles using Equation (5). The right column gives the flow behavior according to the augmented solution. The left bottom image in Figure 17 shows the problem arising when the prior solution is at a large angle with respect to the far-field flow: fluid appears to be deflected in a horizontal direction, which is physically unrealistic. The problem does not occur for the 90 • rotated square [bottom right; using Equation (12)], which shows the same flow behavior as for the topmost right image, henceforth the solution for the large incidence angle may be assumed to be correct.

Discussion
The areal doublet is an analytical element formed by superposing an infinite number of line doublets. The mathematical function can be used to model the physical flow of fluids in highly conductive, narrow flow-channels such as natural fractures [5,6]. The areal doublet locally accelerates the flow of fluids due to the enhanced permeability contrast with the matrix and can be superposed with other flow elements such as nearby hydraulic fractures to model the fluid withdrawal patterns (drained rock volume-DRV) in hydraulically fractured hydrocarbon and geothermal reservoirs [8,23]. The work presented here provides a modification of the algorithm for an areal doublet, which corrects the high deviation of streamlines, when the areal doublet is not parallel to the direction of the fluid flow. The original equation (Equation (5)) results in the same solution as the augmented algorithm (Equation (12)) when the flow direction is parallel to the natural fractures as shown by Figures 9a and  13a. However, the DRV becomes progressively unrealistic when the inclination of the areal doublet increases with respect to the far-field flow (Figure 9b-f). The augmented algorithm superposes two doublet solutions with different vertices based on the tilt angle to generate realistic results for doublets oriented at any angle with respect to the far-field flow. The superposition was needed to avoid the occurrence of branch cuts in inconvenient locations, while still providing corrections to the prior solution. Some key observations from the current study are briefly discussed below.

Effect of Permeability Contrast
The augmented solution corrects the high deviation of streamlines that occurs when the natural fractures are not oriented parallel to the direction of the flow. The augmented solution has a significant impact on results when the permeability contrast is high, as shown by the comparisons of Figures 15  and 16. The permeability contrast for Figure 15c,d is 10 times higher than the permeability contrast of Figure 15a,b ( Table 2). The augmented solution (Figure 16a-d), shows that the DRV reaches slightly further into the reservoir than before (Figure 15a,d), which may increase interference between adjoining wells. The high permeability contrast with the matrix results in larger distortion of the DRV and the impact increases when natural fractures occur at oblique angles with respect to the far-field flow. The DRV is significantly shifted from the original location, whereas the effect of natural fractures that are either parallel or perpendicular to the far-field flow is much more muted. The orientation of the natural fractures, the permeability contrast of the individual natural fractures with the matrix rock, and the natural fracture density, will control the final shape of the DRV.

Model Strengths and Weaknesses
Discrete finite volume methods are widely used to model flow in porous media. Models based on CAM have several strengths, which overcome some of the drawbacks of grid-based methods. The most significant strength of the CAM method is the high resolution, which allows for a closer examination of flow patterns near tightly spaced hydraulic fractures. Another benefit is the high computational speed due to the avoidance of arduous gridding thanks to the closed-form formulation of algorithms. Models based on CAM have been validated with commercial numerical reservoir simulators to verify the accuracy of the results [25,26]. CAM studies allow differentiation between the tracer front (calculated from the outer boundary of the DRV) and pressure front (calculated from the diffusivity equation) [27,28]. Previously, it was assumed that the pressure front marks the drainage boundary around a horizontal well and a hydraulic fracture, but the results from the CAM model showed that pressure front highly exaggerates the actual DRV.
As with any new technology, CAM has some limitations and drawbacks. Current modeling is only limited to 2D single-phase flow. Although, some efforts have been made to expand the model to 3D [23], the additional complexity results in loss of computational efficiency.

Conclusions
A modification to an algorithm developed in earlier studies [5,6] is presented to account for the unusual deviation of the streamlines when natural fractures (represented by the areal doublets) are oriented at large angles with respect to the direction of a far-field flow. The algorithm presented earlier [5] had been validated by using an independent numerical simulator when the flow direction is parallel to the areal doublets or at a moderate angle [25,26].
Our earlier paper [6] was concerned with the evasion of branch cuts such that the flow across natural fractures could still be modeled by continuous Eulerian particle tracking, avoiding discontinuities by smart placement of the branch cuts. In the present paper, the solution of the previous paper [6] was augmented in several steps. We first derived Equation (8), which can give the correct instantaneous streamline solutions by the integral method, but not by the Eulerian particle tracking method. Solutions for the Eulerian method are needed because it allows (1) particle path solutions for both steady and time-dependent flows, and (2) time-of-flight contouring for fluid withdrawal studies of hydrocarbon resources. Subsequently, several additional steps resulted in the augmented solution of Equation (12), which sidesteps the branch cut issues and can be used with both the integral method for instantaneous streamline solutions and with the Eulerian particle tracking method. The augmented solution is an important achievement, because it can accurately model the flow paths across discrete natural fractures, regardless of their orientation with respect to the far-field flow. The prior solution in the earlier branch cut solution paper [6] would become progressively inaccurate for high-incidence angles of the far-field flow to the natural fractures.
The current augmented solution was validated by rotating a rectangular region into different alignment directions with the far-field flow direction (Section 4.5). Based on the results in the current study, the following conclusions can be drawn: a.
Areal dipole/doublet solutions suffer from branch cuts due to the multivalued nature of the algorithm. The branch cuts can be manually removed, but the procedure is cumbersome [6,10]. b. Areal doublets can be used to model the flow in high conductivity flow channels such as the natural fractures in porous media. The natural fractures (areal doublets) distort the DRV and locally increase the velocity while modeling the flow in porous media (Figure 16). c.
The modified areal doublet algorithm in Equation (12) yields a more realistic result for fractures oriented at a high angle with respect to the direction of fluid flow (Figure 16d). d. The augmented solution in the present study avoids the branch cuts and therefore provides an efficient method to model the flow paths of fluids in naturally fractured porous media.