Modeling Flow and Pressure Fields in Porous Media with High Conductivity Flow Channels and Smart Placement of Branch Cuts for Variant and Invariant Complex Potentials

A long overdue distinction between so-called variant and invariant complex potentials is proposed here for the first time. Invariant complex potentials describe physical flows where a switch of the real and imaginary parts of the function will still describe the same type of physical flow (but only rotated by π/2). Such invariants can be formulated with Euler’s formula to depict the same flow for any arbitrary orientation with respect to the coordinate system used. In contrast, variant complex potentials, when swapping their real and imaginary parts, will result in two fundamentally different physical flows. Next, we show that the contour integrals of the real and imaginary part of simple variant and invariant complex potentials generally do not generate any discernable branch cut problems. However, complex potentials due to the multiple superpositions of simple flows, even when invariant, may involve many options for selecting the branch cut locations. Examples of such branch cut choices are given for so-called areal doublets and areal dipoles, which are powerful tools to describe the streamlines and pressure fields for flow in porous media with enhanced permeability flow channels. After a discussion of the branch cut solutions, applications to a series of synthetic and field examples with enhanced permeability flow channels are given with examples of the streamline and pressure field solutions.


Introduction
Complex potentials have been used extensively to describe physical flows [1,2]. The complex potential comprises the potential function as the real term and stream function as the imaginary term, which can be used to describe the pressure field and the velocity field, respectively. The potential function and the stream function satisfy the Laplace equation; thus, the contours of the solutions of the velocity potential and stream function are perpendicular everywhere in the flow space. Switching of the imaginary and real parts of certain complex potentials provides solutions for physically different flows, while the same switch operation in other complex potentials does not change the flow physics. The complex potentials of distinctive flows are coined here variant, as opposed to invariant complex potentials where the switch of imaginary and real axes results in complex potentials that describe the same flow pattern, but only rotated by π/2.
The distinction between variant and invariant complex potentials introduced in the first part of our study is particularly useful when analyzing flow descriptions created by the superposition of various complex potentials. For example, a doublet flow results from the superposition of a source and a sink flow-The latter each described by variant complex potentials-whereas the doublet's complex potential is invariant (see later). Multiple doublets were superposed in a recent study [3] to construct a new analytic element, suitable to model flow in fractured porous media based on complex analysis methods. The algorithm for the analytical element of flow in fractures enabled fundamental models [4,5] and applied models [6][7][8][9] that aim to advance our understanding of flow in fractured porous media.
The earlier treatise by Van Harmelen and Weijermars [3] left undiscussed (and unresolved) the branch cut choices of the potential function, which are relevant for pressure field solutions. The present study offers an extension of the earlier study [3] and allows calculation of the pressure field after considering various branch cut solutions. The solution path presented here is practical for use in applications that require quantification of both the pressure field and the particle paths in a fractured porous medium. In this study, we also present several other basic examples of complex potentials that are either variant or invariant, depending on whether the switching of the imaginary and real parts of the complex potential results in physically distinct flows (variant complex potentials) or the same flows (invariant complex potentials).
This work proceeds as follows. Section 2 first explains the fundamental difference between invariant and variant complex potentials, starting with simple flows (Section 2.1 and 2.2) and then expanding to complex potentials of superposed flows (Sections 2.3 and 2.4). Section 3 explains how the earlier solution using multiple doublets to model flow in fractured porous media [3] is an invariant complex potential. Various options are given to place the branch cuts in the most convenient locations when generating pressure field solutions. Section 4 provides an application to fractured porous media by calculating the stream function and potential function solutions for synthetic and natural examples of fracture sets. Section 5 is a discussion, and Section 6 gives conclusions.

Methodology: Swapping Complex Potentials
The complex potential, which has the potential function as the real part [ ( , )] z t ℜ Ω The real and imaginary parts of certain single-valued complex potentials provide physical descriptions for respectively the potential field contours or isobars (see later Equation (10a)) and the flux across specific regions in the flow between specific streamlines. Although complex potentials for 3D flows are possible [10][11][12], the most convenient formulations are for 2D flows.
A distinction is proposed here for the first time, between so-called variant and invariant complex potentials. Invariant complex potentials describe physical flows where a switch of the real and imaginary parts of the function will still describe the same flow (but only rotated by π/2). Such invariants can be formulated with Euler's formula to depict the same flow for any arbitrary orientation with respect to the coordinate system used (see Section 5.1). In contrast, variant complex potentials, when swapping their real and imaginary parts, will result in two fundamentally different physical flows. Table 1 lists examples of invariant and variant complex potentials, which are illustrated in the present study. Although the areal dipole is categorized as an invariant complex potential in Table 1, a more nuanced study later (Section 3.3) shows that branch cuts may make them variant. All flows studied occur in a 2D flow domain with unconfined boundaries.

Simple Invariant Complex Potential: The Far-Field Flow
The complex potential for a uniform flow or a far-field flow of velocity parallel to the x-axis, x u (m 2 ·s −1 ), is: where z = x + i·y. The potential function (φ ) and the stream function (ψ ) at every point of the flow space can be mapped by the real and the imaginary parts of the complex potential ( ) φ ψ Ω = + z i , as shown in Figure 1a,b. When the real and the imaginary parts of the complex potential of Equation (2a) for the far-field flow are switched, the new complex potential becomes: For oblique flows, one may use i e θ instead of I for the orientation factor in Equation (2b).
The potential function (φ ) and the stream function (ψ ) in every point are mapped by solutions for the real and the imaginary parts (Figure 2a,b) show that Equation (2b) still describes the same flow as in Equation (2a), except for a rotation of the coordinate system, but rotated by π/2. Consequently, the complex potential for a far-field flow is termed here invariant (Table 1).

Simple Variant Complex Potential: The Source/Sink and Vortex Flows
This section shows how source/sink flows (Section 2.2.1) and a vortex (Section 2.2.2) are described by similar complex potentials, but nonetheless are physically different flows; their complex potentials are variant. A final subsection (Section 2.2.3) introduces the final choices of branch cuts to be made, even in the case of a simple variant complex potential (the vortex).

Source/Sink Flows
The complex potential for a source/sink flow of strength, m (m 2 ·s −1 ), centered about an arbitrary location, zd is: The potential function, φ , in every point is solved by the real part of the complex potential ( ) z i φ ψ Ω = + . Switching to polar coordinates (r,θ ), and using identity ln ln z r iθ = + , gives the potential function ( Figure 3a): The stream function, ψ , is given by the imaginary part of the complex potential ( Figure 3b): It is well known that the potential function and stream function for a vortex are described by inverting the real and imaginary parts of the complex potential for a source/sink flow [1]. The latter can be obtained from the complex potential for a source/sink flow through the division of the righthand term of Equation (3a) by an imaginary number i, and replacing the source strength, m (m 2 ·s −1 ), with Γ (m 2 ·s −1 ), denoting the fluid circulation around the vortex singularity [1]: The vortex circulation is defined by vdr Γ =    , such that a positive Γ gives counterclockwise flow and negative Γ results in clockwise flow. The pressure field potential, φ , in every point is solved by the real part of the complex potential, again using polar coordinates (r, θ ) for simplicity ( Figure 4a): The stream function, ψ , given by the imaginary part of the complex potential ( Figure 4b):  Figure 2a). However, the complex potentials of Equations (3a) and (4a) clearly describe two flows that are physically different. Such complex potentials are termed here variant (Table 1).  Flow is anti-clockwise from low to high potentials; pressure gradient will be opposite to the potential gradient.

Branch Cuts
In addition to the distinction of variant and invariant complex potentials, our study also pays particular attention to branch cut choices. Branch cut discontinuities occur when either the potential function (real part of a complex analytic function) or stream function (imaginary part of a complex analytic function) is multi-valued. Usually, multi-valued solutions do not occur simultaneously for both the potential and stream functions of a given flow. The discontinuity appears as a branch cut when the function is visualized in the complex plane, which accommodates the value jump. The branch cut makes the valid solution domain simply connected and the function is then single-valued. A modified contouring, which involves concatenation of different branches is also presented in Section 5.2. An alternative way for representing multi-valued functions uses Riemann surfaces instead of branch cuts.
The stream function solution of Figure 3b is not a branch cut that would pose any challenge in physical interpretation, because the segmented region where the minimum and maximum of the stream function meet does not pose any physical barrier for the flow, only an arbitrary reference line for the flux computation. However, the potential function solution of Figure 4a we label as a critical branch cut choice, which would pose an apparent obstacle for 2D flow across the branch cut. Fluid in the clockwise rotating vortex flow would need to jump over the implied pressure gradient barrier at the location of the branch cut, from low to high pressure. This is purely a mathematical artefact or choice that assigns a single value to a multi-valued complex elementary function [13], which nonetheless could create confusing plots if not properly interpreted. For the vortex field, any radial may be chosen as the desired branch cut location (Figure 5a-d), which highlights that the vortex has a quasi-continuous potential function (pressure gradient). However, when we analyze flows due to complex potentials of superposed flow elements, the location of the branch cuts may leave open many choices, some less practical than others, as we will discuss in detail in Section 3.3.

Superposed Invariant Complex Potentials: Doublet Flow
The complex potential for a so-called spaced doublet is derived by combining the vector fields for a point source and a point sink of equal strength located at z1 and z2 respectively (Figure 6a). For a singularity doublet, a limiting process that decreases the distance between the point source and sink (z1-z2) to zero (Figure 6b), gives the required complex potential [14]. The singularity doublet and singularity dipole can both be described by an invariant complex potential (Sections 2.3.1 and 2.3.2). The distinction between doublets and dipoles was originally proposed by Strack [15], which was required for constructing line arrays called, respectively, line doublets and line dipoles. The distinction is made based on the polarity orientation of these elements (for later use in a superposition procedure. see Section 2.4). The line doublet has the unique straight streamline of each singularity doublet oriented perpendicular to the line array. In contrast, the line dipole aligns the symmetry streamline with the line array of all singularities. The distinction between doublets and dipoles by Strack [15] was adopted in a later study by Weijermars and van Harmelen [14], and is further analyzed in the present study.

Singularity Doublet
The complex potential for a singularity doublet with a strength of µ (m 3 ·s −1 ) located at zd in a complex plane, obtained by integrating the vector field is [14]: .
The pressure field potential, φ , is the real part of the complex potential ( ) for a singularity doublet, using polar coordinates (r,θ ) for simplicity (Figure 7a), is: The stream function, ψ , due to the imaginary part of the complex potential (Figure 7b), is:

Singularity Dipole
When the real and the imaginary parts of the complex potential of Equation (5a) for the singularity doublet are switched, through division of the right-hand terms in Equation (5a) by an imaginary number i, and replacing the doublet strength, µ (m 2 ·s −1 ), with µ* (m 2 ·s −1 ), the new complex potential becomes: The potential function, φ , and the stream function, ψ , for the rotated point doublet ( Figure   8a,b), using polar coordinates (r,θ ), are: The complex potentials of Equations (5a) and (6a), if scaled by the same strength µ = µ*, clearly describe two flows that are physically identical (save a rotation over π/2). Such complex potentials are termed here invariant (Table 1). Although the rotated singularity doublet flow of Figure 6b is physically identical to that of Figure 6c (except for a reference frame rotation), the flow of Figure 6c is termed a point dipole, to aid the distinct of the superposition series in Section 2.4.

Superposed Variant Complex Potential: Line Doublets and Line Dipoles
A point doublet, or singularity doublet, is an analytical element where one side of the singularity is the source (+) and the other side is the sink (−) (Figure 6b). The point (or singularity) dipole is an analytical element like the singularity doublet, but the source and sink sides are rotated 90° clockwise ( Figure 6c). The distinction between a singularity doublet and singularity dipole becomes important when considering a linear array of superposed singularity doublets ( Figure 9a) or a linear array of superposed dipoles (Figure 9b). The flow regimes of each of these superposed arrays are physically different and have distinct implications in hydrocarbon well placement, geothermal reservoirs, and flow through a heterogeneous reservoir [3,14,15].

Line Doublets
The complex potential for modeling the flow near discrete fractures makes use of the superposition of line doublets (see Section 4), and can be obtained from the complex potential for a singularity doublet with a strength of m (m 3 ·s −1 )located between za and zb in a limiting superposition process (Figure 9a): The potential function (Figure 10a), φ , is mapped by the real part and the stream function ( Figure 10b), ψ , is given by the imaginary part the complex potential ( ) respectively.

Line Dipoles
The complex potential for a line dipole with a strength of m (m 3 ·s −1 ) located between za and zb in a complex plane, can be obtained by superposing dipole singularities in a limiting process ( Figure  9b): The potential function (Figure 11a), φ , and the stream function ( Figure 11b), ψ , are mapped by the real part and the imaginary part of the complex potential respectively (see Equation (1)). The complex potentials of Equations (7) and (8) clearly describe two flows that are physically different. Such complex potentials are termed here variant (Table 1).

Complex Potentials to Describe Flow in Enhanced Permeability Flow Channels-Invariant Complex Potentials
Van Harmelen and Weijermars [3] introduced an analytical element to describe flow through permeable fractures, which enhance fluid velocity locally, by superposing an infinite number of line doublets. Figure 12a shows the spatial arrangement of line doublets to be integrated to obtain what is coined here as an areal doublet element. The same element was originally referred to as a natural fracture element [3], because the element may model flow in discrete fractures. However, using the term "areal doublet element" in hindsight seems more objective than "natural fracture element" as the element can model flow accelerators in porous media. Below we show that an element similar to the areal doublet can be obtained via an alternative superposition of line dipoles, following the spatial integration method sketched in Figure 12b.

Branch Cut Issues
Before discussing the areal doublet, Figure 13a What is potentially confusing is the occurrence of an "apparent" adverse pressure gradient when the number of superposed line doublets increases, and the branch cut spacing is reduced accordingly. The pressure gradient in the central region between the outer line doublets remains favorable to flow from high to low (as opposed to adverse pressure gradients), but involves many pressure jumps. When the number of branch cuts increases, a superficial inspection may suggest and adverse pressure gradient develops in the central region between the outer branch cuts. This flow adversity of the potential function gradient is purely a resolution effect, the pressure gradient remains favorable to flow, but the difference between branch cut jumps and pressure gradients is no longer visually discernable when the number of superposed line doublets increases.
The branch cuts are "folded" into the fracture zone confined between the vertices za1-za2 and zb2-zb1, using the vertices as defined in Figure 14. In fact, there is an infinite number of branches due to the analytical element being formed by integrating numerous discrete doublet elements. The inner branch cuts get cancelled as seen in Figure 13a-d, and only two outer branches remain (at the entry and exit of the flow channel). This is similar to what is observed in Figure 13a-d where only the other two branch cuts remain in the potential field when the number of line doublets are increased.

Areal Doublet
The complex potential for the areal doublet, used in a prior study of flow in fractures [3], is obtained by superposing an infinite number of line doublets (see Equation (B12) in [3]): where υ (m 4 ·s −1 ) is the strength, which is assumed to be steady here, but may be time-dependent for certain cases, γ (radians) is the tilt angle with x-axis, β (radians) is the angle between the corner points ( Figure 14), h (m) is the height, L (m) is the length, W (m) is the width, za1, za2, zb1, and zb2 ( Figure 14) are co-ordinates for the vertices calculated from center co-ordinate, zc, as follows: The tilt angle with x-axis, β (radians), in this study is exclusively used as π/2 in order to create the flow channels which are rectangular in shape. The pressure function (Figure 15a,b), φ , and the stream function (Figure 15c), ψ , are mapped by the real part and the imaginary part of the complex potential respectively (see Equation (1)).

Branch Cut Choices for Areal Doublet
The solutions of the potential function (Equation (9c)) and stream function (Equation (9d)) by standard software platform solvers introduce branch cut choices. For simpler expressions (such as solutions of Figure 5a-d), the solver functions embedded in the mathematical software work well to separate the real and the imaginary parts of the complex potential to plot the potential and the stream functions, respectively. However, for an expression as complex as the one in Equation (9a), the expression must be manually solved to avoid the appearance of branch-cuts in impractical locations (see Appendix A, Appendix A3). Several logarithmic functions, with non-zero argument, are present in Equation (9a)  The pressure fields of Figure 15a,b are meant to illustrate the problem that arises with the scaling of the outer pressure field (Figure 15a, when not folding the branch cuts into the fracture) due to the pressures magnitude monotonically increasing to an infinitely large value in the zone between the branch cuts which extends to infinite (outside the top of the picture). By folding the branch cuts inside the fracture (Figure 15b), the pressures values immediately outside the fracture will not be suppressed by infinitely high pressures arising in Figure 15a as a result of the unfavorable branch cut placements.
Contour integration of several mathematical functions involving square roots, power, natural logs, and inverse trigonometric functions, are multi-valued and thus may result in branch cuts appearing in inconvenient locations when plotting integral solutions. The branch cuts may occur in pressure field solutions as well as in stream function solutions for various closed-form solutions, as seen in Figures 3b, 4a, 10a and 11b. These branch cuts may mask the intended values of the physical quantity (such as pressure) described by the complex potential. The location of the branch cuts can be better controlled when the potential function is manually separated into the real and imaginary parts, using Equation (9a) and following the schedule explained in Appendix A (Appendix A3). Figure 15b shows a solution where the branch cuts are not extending toward infinity, but are pointing inward, confined to the central area (Figure 15b), similar to the solution obtained in Figure 13d (after adding an infinite number of line doublets in the superposition). The solution of Figure 15b suggests the existence of an adverse pressure gradient in the central flow region, which is entirely due to the method of solution, similar to what is seen in Figure 13a

Areal Dipole
The complex potential for superposed line dipoles can be obtained from the complex potential for a superposed line doublet (Equation (9a)) through the multiplication of the right-hand term of Equation (9a) by an imaginary number i. The pressure field potential and stream function is mapped by the real and imaginary part of the complex potential respectively as shown in Figure 16a,b. The branch cut appearing in the stream function solution of Figure 16b is moved using the same approach as for the branch cut relocation of the potential function in Figure 15a. The conclusion is that the complex potentials for both the areal doublet and areal dipole are invariant (Table 1), because switching of its imaginary and real parts (or equivalent swapping of the real and imaginary coordinate axes) gives their mutual solutions. The complex potentials for the areal doublet and areal dipole describe two flows that are physically identical. What is remarkable is that the line doublets and line dipoles that are used to obtain the invariant complex potentials of the areal doublet and areal dipole, are themselves variant (Table 1). However, the selection of branch cuts, based on the nature of the problem, may result in areal doublet to be variant.

Application
This section is sub-divided into several subsections. The basic tools for pressure scaling, streamline tracking and time of flight contours are outlined in Sections 4.1 and 4.2, respectively. Sections 4.3 and 4.4 apply the developed tools to model several synthetic flow elements, including high conductivity flow channels with randomized orientations, lengths, and strengths, which may mimic natural fractures. Finally, Section 4.5 analyzes the flow of fluid through a non-isotropic system with several high conductivity flow channels abstracted from a field example.

Scaling of Pressure Field Solutions
The transient pressure gradient, P(z,t), governing the time-dependent fluid flow in a reservoir with permeability k can be calculated from the real part of a suitable complex potential, The complex potential is split in its real and imaginary parts, which correspond to respectively the potential function, ( , ) z t φ , and the stream function, ( , ) z t ψ as shown in Equation (1). Subsequently, the potential function, ( , ) z t φ , can be used to compute and scale the pressure gradient due to a flow (relative to an initial reservoir pressure) in any location, z, at a given time, t, using the following generic formula [16,17]: where ( , ) z t φ is the potential function, η the viscosity of the fluid (expressed in Pa.s), and k the permeability of the reservoir (expressed in m 2 nor Darcy). The generic form of Equation (10b) will yield a specific pressure profile when a specific potential function is defined (and η and k are known).
The actual pressure field at any given time, t, can be thus be obtained when the initial pressure state, P0, of the reservoir is known: The porosity may appear as a scaling factor accounted for in input values for ( , ) z t φ such as a certain flow strength [18].

Particle Tracking and Time-of-Flight Contours
The algorithm used to trace streamline and time-of-flight contours apply a Eulerian particle tracking method ( Equation (11) is used to calculate the velocity field solutions from specific velocity field expressions defined for any analytical flow elements. Tracing each streamline is accomplished by first choosing an initial position z0 at time t0 = 0 and calculating the initial velocity. By choosing an appropriate time-step, Δt, the position of the tracer at time t1 = t0 + Δt is: The position of the tracer at any other time tj is: The time-of-flight contours (TOFCs) are determined by plotting the location of all the tracers for particular time steps. The selection of an appropriate time-step duration is a crucial decision in the Eulerian particle tracking method. An overly coarse time step in combination with sharply curving streamlines results in overshooting of the particles, which then jump from one streamline to the other, which is inaccurate. For such cases, a very small time-step is required such as to avoid inaccuracies in the displacement of the particles due to the high instantaneous velocity. However, for cases where the velocity gradient is spatially less convolute, a coarser time-step can be used to reduce the computational time. Further details on complex potentials and stream functions in reservoir flow applications are given elsewhere [2][3][4][5][6][7][8][9]16,17].

Synthetic Cases of Wide Flow Channels (Cases A-B)
A first series of flow channels is modelled with large apertures or widths to highlight the streamline convergence in the flow element as well as the fact that the steepest pressure gradients occur in and near the flow elements. For each of the simulations in this section, a uniform far-field flow rate of 0.15 m·year −1 (4.75 × 10 −8 m·s −1 ) is assumed, from left to right parallel to the real axis. For later comparison, Figure 17a first shows the streamlines (blue) and TOFC (red), due to the far-field flow in an isotropic reservoir without any other flow elements (such as hydraulic/natural fractures) simulated for the period of 20 years. The fluid viscosity is assumed to be 1 cP (0.001 Pa.s) and the reservoir permeability is 10 mD (9.87 × 10 −15 m 2 ). The particle tracking simulation is run for a period of 20 years with a time step (Δt) of 0.1 days. For the far-field flow travelling at a rate of 0.15 m·year −1 , the streamlines and the TOFC travel 30 m downstream (in the homogeneous and non-fractured reservoir). The porosity is assumed to be 10%; for cases with different porosity the TOFC will be proportionally slower or faster (see [18] for additional details).   Table A4). (c) Case B: Two enhanced permeability flow channels perpendicular to each other (attributes given in Appendix B, Table A5). Right column: Resulting pressure fields in Pascal.

Effect of Spatial Complexity in the Location of the Flow Channels (Cases C-E)
The next series of simulations use a narrower and increasing number of enhanced permeability flow channels, which highlights that when the number of flow channels increases (1) the pressure gradients become steeper, and (2) the pressure contour pattern becomes more complex. For the base case, consider a 60 m × 60 m unconfined isotropic flow space (Figure 18a). A uniform far-field flow moves from left to right with a velocity of 0.15 m·year −1 (or 4.7 × 10 −8 m·s −1 ), which maintains a maximum pressure of 2.9 × 10 2 kPa for a reservoir with a porosity of 10%. Figure 18b-d show the pressure fields due to the presence of increasing numbers of flow channels. The flow channels are randomly generated as follows. The width of each flow channels is randomly varied between 0.1-0.01 m, a length between 5 and 20 m, and orientation between −90° to 90° with respect to the horizontal axis. The strength of the is varied between 4.7 × 10 −8 and 9.4 × 10 −6 (m 4 ·s −1 ). The maximum pressure and the local pressure gradients both increase when the number of flow channel increases. Simultaneously, the pressure contour pattern becomes more complex.

Field Example from Aravalli Supergroup-Rajasthan (India) (Case F)
A final set of simulations revisits flow along the hydrothermal veins [3] using the fracture pattern observed in a polished slab of facing stone quarried from fractured Proterozoic rocks from the Aravalli Supergroup in the state of Rajasthan, India [19][20][21]. The fractures are formed under high fluid pressure deeper in the Earth's crust before the rock was exhumed by tectonic uplift and erosion. These rocks are currently exploited as facing stones and quarried near the villages of Bidasar-Charwas, Churu district (Figure 19a). The quarries are confined to a 0.5 km wide and 2.5-3.5 km long belt of open pits dug below the desert plain. The rocks in these pits have been described as part of the Bidasar ophiolite suite [22]. A polished slab containing the naturally created fracture network are imaged in Figure 19b. Similar to our companion study [3], the simulations in this study do not account for the creation of the fractures, which are already developed when invaded by the hydrothermal fluids. The precise natural pressure responsible for the injection of the hydraulic veins is unknown, but the pressure has exceeded the strength of the rock and was large enough to inject fluid into propagating fractures at several km burial depth, thus being in the order of 100 MPa. Figure 19c highlights the main flow conduits used in our flow analysis. The far-field flow velocity is assumed to be 4.75 × 10 −9 m s −1 and porosity is 10% for all rock. The direction of the far-field flow in Figure 20a-d is assumed to be from top to bottom (−90°), based on the provenance of the major flow channels. Figure 20a shows a base case for comparison with the time-of-flight contours and streamlines in the absence of the enhanced permeability flow channels.   Table A8).
The pressure field plots for each case (Figure 21a-d) shows that, for cases with higher permeability contrasts (Figure 21c,d), the maximum pressures increase and so do the local pressure gradients, while the pressure contour pattern becomes more involved. Figure 21c,d also show a highpressure zone develops in the left part of the flow field where numerous significantly long horizontal and vertical flow channels are aggregated.
The pressure field generated in Figure 21c,d show anomalously high and low pressure zones in some areas of the reservoir (which can be termed as hot/cool spots). These hot-spots and cool-spots are the artefacts of the arbitrarily specified strengths of the flow channels and occur for the cases where the natural fracture strengths are unusually high (Figure 21c,d). The effect does not appear for cases where the natural fracture strength is reasonable (Figures 18b,c and 21b). The pressure near the flow channels are also distorted as the potential function is scaled by a uniform permeability (Equation (10a)), while in the model the flow is artificially enhanced by the analytical element. The pressure field needs to be downscaled in the vicinity of the flow channel elements by a locally decreasing the permeability. The realistic scaling of the pressures in each flow channel and selective downscaling of the permeability inthe vicinity of the fractures requires further study (future work).

Variant and Invariant Complex Potentials
The work presented here differentiates between two types of complex potentials, the variant and the invariant complex potential. For several complex potentials, swapping of the real and the imaginary axes (or the equivalent operation of multiplying the potential by an imaginary number i), results in a flow which is physically different from the flow represented by the original complex potential. Examples of such variant complex potentials include the point source ( Figure 3) and vortex (Figure 4), line doublet ( Figure 10) and line dipole (Figure 11), which are fundamentally different flows. For other flows, the complex potential remains invariant, because even when axes are swapped, the same flow occurs, only rotated by π/2. Examples of invariant complex potentials include the far-field flow (Figures 1 and 2) and the singularity doublet (Figures 7 and 8). For the farfield flow the invariance of flows in any direction is apparent from: The invariant nature of complex potentials can for simple flows be exploited to obtain the stream function solution from the potential function, and vice versa. However, for complex flows with invariant complex potentials, further operations may be required to compute the stream function or potential function. Table 1 categorizes the areal dipole and areal doublet as invariant complex potentials. However, based on further discussion presented in Sections 2.2.3, 3.1 and 3.3 (regarding branch cuts), only certain branch cut choices for the areal dipole/doublet show the invariant flow outcome. For example, solutions with branch cuts obtained by separating the real and the imaginary parts manually, as shown in Figures 15b and 16c, are invariant. Other placement of branch cuts (see later in Section 5.2) would render the result variant.

Trade-Offs in Branch Cut Choices
The superposition of an increasing, discrete number of line doublets in Section 3.1 to create the areal doublet (Figure 13a-d) showed that the number of branch cuts increases commensurate with the number of superposed line doublets. The simplest branch cut positioning is a straight-line coinciding with each line doublet. When the process is continued to infinity, the entire central area confined between the two outer line doublets will be occupied by branch cuts, meaning that the potential function is not single-valued in the region portraying the fracture. This of course is a mathematical artefact, and the solution of the stream function shows that flow does occur in that region. We undertook various efforts to either repair the branch cut, or move the branch cut due to the potential function to another region. Below we discuss the various options with their limitations and trade-offs.
The potential function solution of Figure 15b suggests a continuous pressure gradient occurs in the central region. This is an arbitrary solution for the region where the potential function is multivalued, suggesting continuity of both the potential lines and gradients across the central areal doublet. However, we emphasize this smoothed graphical solution is physically contradictory with the flow direction through the central area, which is counter to the pressure gradient. The stream function solution for the areal doublet does not involve any branch cut and is everywhere singlevalued and therefore physically correct. Although adverse pressure gradients may arise near boundary layers in flows involving inertia effects [23], no such effects can occur in the Darcy flows through the porous media considered in our study. Our sub-conclusion regarding the potential function solution of Figure 15b is that the central region is completely occupied by branch cuts, and the continuity of the contour gradients is an arbitrary solution, counter-intuitive because of the creation of an apparent adverse pressure gradient, which is non-physical with regard to the stream function solution (Figure 16c). Figure 15a shows an alternative solution, where the branch cuts for the potential function occur exterior to the central areal doublet region. This solution shows a realistic potential gradient inside the fracture element but generates branch cuts in the exterior region that extend to infinity. This solution with exterior branch cuts becomes unworkable and impractical when superposing numerous areal doublets in the far-field flow to simulate multiple fractures (Section 4). Exterior branch would obliterate part of the pressure field around the fractures. Placing branch cuts inside the fractures is more practical (Figure 15b), especially when these are narrow in width, because these occupy a limited area relative to the matrix. Below we explain other examples of moving the exterior branch cuts to the interior of the areal doublet. Figure 22a shows another alternative solution, where the branch cut is placed at one end of the areal doublet region. This give a physically plausible and so-called favorable pressure gradient for the flow region in the left half of the image where no branch cut occurs, but still is impractical for the right half of the image. The final solution we offer is a composite solution, where the left and right half of the flow region is free from branch cuts and with just one branch cut occurring centrally in the doublet flow region (Figure 22b) as explained in Appendix A (Appendix A3, last paragraph) which is similar to sequential contouring presented by Holzbecher [24]. This solution also suppresses the extreme highs and lows in the pressure scale that occur when the branch cuts are pointing outward to infinity (Figures 15a and 22a). Although the potential function solution of Figure 22b is physically the best choice, it is not practical in execution. For pressure plots with multiple fractures we use the solution of Figure 15b, which is computed after separation of the real and imaginary parts of the complex potential as explained in Appendix A (Appendix A3).

Applicability and Accuracy of Current Approach
This study uses an alternative method to the widely used technique of numerical simulation with discrete volumes to model the flow of fluid in porous media with highly conductive flow channels. Unlike such discrete volume models, the present method based on complex analysis methods is grid-less, mesh-less and much faster than the discretized models. Additionally, the closedform solution provides infinite resolution. The Eulerian particle tracking method implemented in this study gives unique insights by visualizing the flow in fractured reservoirs at high resolution ( Figures  17 and 20).
The pressure field plots for the areal doublet used to model the flow and pressure field near and in the enhanced permeability flow channels, involve choices about the placement of branch cuts to curtail multi-valued solutions, which will locally affect the minimum and maximum pressures. However, the pressure field outside the enhanced permeability flow channels show consistent singlevalued results (Figures 17, 18, and 20). The width of these flow channels is generally in the range of micro-to millimeters [25,26], so for practical purposes, the pressure estimations presented in this study will introduce only minor inaccuracy due to the branch cut choices made. Although we don't assert that the areal doublet is a mathematical proxy for natural fractures, they can certainly model flow elements which locally enhance fluid velocity during a flow in porous media.

Conclusions
In this work, we present several analytical flow elements which can be used to represent the physical flows. For some of these flows, switching of the real and imaginary parts of the complex potential, which can be achieved by multiplying or dividing the complex potential by an imaginary number i, results in a different flow, termed in this study as variant complex potentials. However, for certain other flows, multiplication or division of the complex potential by the imaginary number, i, does not change the flow physics, which are termed in this study as invariant complex potentials. We also present new solutions to calculate the pressure field based on the complex potential for an areal doublet. From the results in this study, we can conclude that areal doublets, which can be used to model high conductivity flow channels, locally increase both the fluid velocity and the pressure gradients (Figures 17 and 20

Conflicts of Interest: Page: 22
The authors declare no conflict of interest.

Appendix A. Branch Cuts in Pressure Plots
For areal doublets (superposed line doublets), simply taking the real part of the complex potential may lead to inconveniently placed branch cuts, which cause discontinuity in the pressure plots. Thus, the real and complex potential need to be manually separated in order to facilitate choices about branch cut placement. Examples for branch cut solutions are presented below.

Appendix A1. Background
The real part of the complex potential yields the potential function, which can be used to calculate the pressure field in a reservoir (as detailed in Section 4.1). Several of our prior studies [12,13] include pressure field solutions for several analytical elements (such as point sources and sinks representing vertical injectors and producers, as well as line sinks, representing hydraulic fractures) and obtained excellent matches with independent pressure field solutions generated using a numerical reservoir simulator.
The superposed complex potential (ΩSUP(z)) for a far-field flow, ΩFF(z) (Equation (2a)), and an areal doublet, ΩAD(z) (Equation (9a)), is: connecting two branch points at ( ) f z = 0 and infinity, which separates the discontinuity present in a complex plane [23]. Next, we calculate the pressure field for a synthetic reservoir with properties as summarized in Table A1 to show the branch cuts for an areal doublet. We assume 0 P defined in Equation (10b) to be zero, which is an acceptable assumption for a simple synthetic case. The corresponding pressure gradient (delta P) can be solved using Equations (A2) and (A3):  Figure A1 shows the pressure profile for a steady-state case, where a horizontal areal doublet with finite width (represented by black dashed lines) is superposed by a far-field flow (from left to right). The fracture affects the pressure field in its vicinity, as can be inferred from the deflected isobars close to the fracture. However, an undesirable computational effect occurs beyond the right end termination of the fracture (represented by solid white lines), where the pressure jump continues for an infinite distance in the horizontal direction toward the right. This effect is due to the occurrence of so-called branch cuts in the solution of the potential function (Equation (A2)), which becomes undefined at the vertices of the fracture. The simple model in Figure A1 has four branch points at the vertices of the fracture, which renders the function in Equation (9a)  z . This results in the pressure profile becoming discontinuous at the branch cuts and the pressure gradient inside and outside of the fracture shows a big jump. The branch in Figure A1 can be better understood if we consider a function ( ) log z , mapped in a complex plane as shown in Figure A2. The value of ( ) log z at A, which is infinitesimally close to and above the positive x-axis, differs from that at B, which is infinitesimally close to A but is below the positive x-axis. The function ( ) log z is discontinuous across the branch cut represented by line connecting zero to infinity. A Riemann surface can be used to represent the multiple values of ω by splitting the z plane into n parallel planes. The discontinuity in Figure A1 can be analyzed by taking two arbitrary points 5 + 10i and 5 − 10i and plotting the pressure profile across the vertical cross-section between those points. This is the region with the branch cuts which extends to infinity, from each vertex. The cross-section of Figure  A3a demonstrates the acute jump of pressure across the branch cut, and as we move infinitesimally close to the branch cut, either from above or below it (−0.5i or 0.5i) as shown in Figure A2. The pressure gradient ( P Δ ) jumps from −0.5 × 10 4 Pa to 1.25 × 10 4 Pa. In effect, each pressure plane lies in two different planes of a Riemann surface as mentioned before. However, Figure A3b shows a pressure profile along −5 + 10i and −5 − 10i ( Figure A1), the pressure gradient is smooth and no discontinuity occur as the branch cuts do not extend to negative infinity.
(a) (b) Figure A3. (a) Pressure profile along 5 + 5i and 5 − 5i in Figure A1. Horizontal scale shows distance in meter and vertical scale is pressure in Pa. (b) Pressure profile along −5 + 5i and −5 − 5i in Figure A1. Horizontal scale shows distance in meter and vertical scale is pressure in Pa.

Appendix A3. Proposed Solution to the Branch Cut Placement
The method adopted here to overcome branch cut effects is to separate the real and imaginary parts manually and calculate the phase angles by using tangent function as shown in Equation (A6) below for all the logarithmic terms in Equation (9a). This results in a smooth pressure profile without any unrealistic pressure jumps as shown in Figure A3a.
The following template is used to manually separate the real and imaginary terms and generate the potential plot solution of Figure 15b: The real and imaginary parts of Equation (A6) can be separated as follows: If the pressure is calculated by using the original formulation (Equation (9a)), the branch cuts will cause apparent jumps (discontinuities0 in the pressure filed that appear to extend in unfavorable zones (za1-zb1, za2-zb2, see Figures 14, 22a and A1). By placing the branches between za1-za2 and zb1-zb2 (Figure 22b), we still have discontinuous solutions known as integral branches, but these branches become confined to the inner region of the flow channels and the pressure field around the fractures at least appear continuous. We recall that this is a pure mathematical manipulation, with the rationale for the chosen solution in Equations (A7)-(A15) being this formulation folds the branch cuts inside the fracture element and achieves the desired continuity of the outer pressure field. Other solutions are possible (Figure 22b).
Equation (A4) with terms A, B, C and D as defined in Equations (A10)-(A15) results in a continuous potential function plot (isobars, Figure A4) for the reservoir defined in Table A1. In Section 2.2.3 ( Figure 5), we demonstrated different branch cut solutions for a simple function involving natural logarithm (vortex). A similar procedure can be applied to combine the branch cuts at a single location between the two ends of an areal doublet element as shown in Figure 22b. The potential function is calculated for each vertex of the areal doublet element shown in Figure 14 by dividing the complex plane into two halves at a particular location. Each opposing vortex pair (za1 vs. zb1 and za2 vs. zb2) is assigned opposite branch location, which results in a single branch cut in the previously selected location. Finally, the solutions are spliced together to obtain the solution in Figure  22b.  Figure A1 has been removed by the applied procedure. Length scale in feet and pressure (rainbow colors) in psi.