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

Abstract

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.

1. 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 Section 2.2) and then expanding to complex potentials of superposed flows (Section 2.3 and Section 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.

2. Methodology: Swapping Complex Potentials

The complex potential, which has the potential function as the real part $\mathsf{\Re }[\mathsf{\Omega }(z,t)]$ and the stream function as the imaginary term $\mathsf{\Im }[\mathsf{\Omega }(z,t)]$, can be used to model two-dimensional flows of an incompressible fluid through a porous medium. The com−plex potential is related to the potential function, $\varphi (z,t)$, and the stream function, $\psi (z,t)$ as follows:

$$\mathsf{\Omega }(z,t)=\mathsf{\Re }[\mathsf{\Omega }(z,t)]+\mathsf{\Im }[\mathsf{\Omega }(z,t)]=\varphi (z,t)+i\psi (z,t)({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(1)

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. Invariant and variant complex potentials.

Invariant	Variant
Far-field flow	Source/sink flow
Singularity dipole/doublet	Vortex
Areal doublet	Line doublet
Areal dipole	Line dipole

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.

2.1. 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, $u_x$ (m²·s⁻¹), is:

$$\mathsf{\Omega }(z)=u_{x}z({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(2a)

where z = x + i·y. The potential function ($\varphi$) and the stream function ($\psi$) at every point of the flow space can be mapped by the real and the imaginary parts of the complex potential $\mathsf{\Omega }(z)=\varphi +i\psi$, 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:

$$\mathsf{\Omega }(z)=i.u_{x}z({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(2b)

For oblique flows, one may use $e^{i\theta }$ instead of I for the orientation factor in Equation (2b).

The potential function ($\varphi$) and the stream function ($\psi$) 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).

2.2. 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).

2.2.1. Source/Sink Flows

The complex potential for a source/sink flow of strength, m (m²·s⁻¹), centered about an arbitrary location, z_d is:

$$\mathsf{\Omega }(z)=(m/2\pi )\ln (z-z_d)({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(3a)

The potential function, $\varphi$, in every point is solved by the real part of the complex potential $\mathsf{\Omega }(z)=\varphi +i\psi$. Switching to polar coordinates (r, $\theta$), and using identity $\ln z=\ln r+i\theta$, gives the potential function (Figure 3a):

$$\varphi (r,\theta )=(m/2\pi )\ln r({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(3b)

The stream function, $\psi$, is given by the imaginary part of the complex potential (Figure 3b):

$$\psi (r,\theta )=(m/2\pi )\theta ({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(3c)

2.2.2. Vortex

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 right-hand term of Equation (3a) by an imaginary number i, and replacing the source strength, m (m²·s⁻¹), with $\mathsf{\Gamma }$(m²·s⁻¹), denoting the fluid circulation around the vortex singularity [1]:

$$\mathsf{\Omega }(z)=\frac{\mathsf{\Gamma }}{2\pi i}\ln (z-z_d),({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(4a)

The vortex circulation is defined by $\mathsf{\Gamma }={\displaystyle \oint \overrightarrow{v}dr}$, such that a positive $\mathsf{\Gamma }$ gives counterclockwise flow and negative $\mathsf{\Gamma }$ results in clockwise flow. The pressure field potential, $\varphi$, in every point is solved by the real part of the complex potential, again using polar coordinates (r, $\theta$) for simplicity (Figure 4a):

$$\varphi (r,\theta )=\frac{\mathsf{\Gamma }}{2\pi }\theta ({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(4b)

The stream function, $\psi$, given by the imaginary part of the complex potential (Figure 4b):

$$\psi (r,\theta )=-\frac{\mathsf{\Gamma }}{2\pi }\ln r({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(4c)

Comparison of Figure 3a,b and Figure 4a,b reveals that what are streamlines for the source/sink flow (Figure 3a) are identical to pressure contours for the vortex flow (Figure 4b). Similarly, the pressure contours for the source/sink flow (Figure 1b) become streamlines for the vortex flow (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).

2.2.3. 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.

2.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 z₁ and z₂ respectively (Figure 6a). For a singularity doublet, a limiting process that decreases the distance between the point source and sink (z₁–z₂) 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 (Section 2.3.1 and Section 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.

2.3.1. Singularity Doublet

The complex potential for a singularity doublet with a strength of µ (m³·s⁻¹) located at z_d in a complex plane, obtained by integrating the vector field is [14]:

$$\mathsf{\Omega }(z)=\frac{i.\mu }{2\pi (z-z_d)}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(5a)

The pressure field potential, $\varphi$, is the real part of the complex potential $\mathsf{\Omega }(z)=\varphi +i\psi$, which for a singularity doublet, using polar coordinates (r, $\theta$) for simplicity (Figure 7a), is:

$$\varphi (r,\theta )=-\frac{\mu }{2\pi }·\frac{\cos \theta }{r}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(5b)

The stream function, $\psi$, due to the imaginary part of the complex potential (Figure 7b), is:

$$\psi (r,\theta )=\frac{\mu }{2\pi }·\frac{\sin \theta }{r}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(5c)

2.3.2. 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²·s⁻¹), with µ* (m²·s⁻¹), the new complex potential becomes:

$$\mathsf{\Omega }(z)=\frac{\mu *}{2\pi (z-z_d)}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(6a)

The potential function, $\varphi$, and the stream function, $\psi$, for the rotated point doublet (Figure 8a,b), using polar coordinates (r, $\theta$), are:

$$\varphi (r,\theta )=-\frac{\mu *}{2\pi }·\frac{\sin \theta }{r}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(6b)

$$\psi (r,\theta )=\frac{\mu *}{2\pi }·\frac{\cos \theta }{r}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(6c)

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.

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

2.4.1. 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³·s⁻¹)located between z_a and z_b in a limiting superposition process (Figure 9a):

$$\mathsf{\Omega }(z)=\frac{-i.m[Log(z-z_a)-Log(z-z_b)]}{2\pi (z_a-z_b)}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(7)

The potential function (Figure 10a), $\varphi$, is mapped by the real part and the stream function (Figure 10b), $\psi$, is given by the imaginary part the complex potential $\mathsf{\Omega }(z)=\varphi +i\psi$ respectively.

2.4.2. Line Dipoles

The complex potential for a line dipole with a strength of m (m³·s⁻¹) located between z_a and z_b in a complex plane, can be obtained by superposing dipole singularities in a limiting process (Figure 9b):

$$\mathsf{\Omega }(z)=\frac{-m[Log(z-z_a)-Log(z-z_b)]}{2\pi (z_a-z_b)}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(8)

The potential function (Figure 11a), $\varphi$, and the stream function (Figure 11b), $\psi$, 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).

3. 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.

3.1. Branch Cut Issues

Before discussing the areal doublet, Figure 13a–d show the potential function (left) and stream function (right) solutions obtained by discrete superposition of different number of line doublets. The potential function solution (left) exhibits an increasing number of branch cuts (Figure 13a–d), corresponding to the number of superposed line doublets. Discrete jumps in the value of the potential function (a proxy for pressure gradients), due to the branch cuts, are graphed in the central images. 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 z_a1–z_a2 and z_b2–z_b1, 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.

3.2. 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]):

$$\begin{array}{l}\mathsf{\Omega }(z)=\frac{-i\cdot \upsilon (t)}{2\pi hn\cdot L\cdot W}.e^{-i\gamma }[(z+z_{a2})(\log (-e^{-i\gamma }(z-z_{a2}))-(z+z_{a1})(\log (-e^{-i\gamma }(z-z_{a1}))+\\ (z+z_{b1})(\log (-e^{-i\gamma }(z-z_{b1}))-(z+z_{b2})(\log (-e^{-i\gamma }(z-z_{b2}))]\end{array}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(9a)

where υ (m⁴·s⁻¹) 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, z_a1, z_a2, z_b1, and z_b2 (Figure 14) are co-ordinates for the vertices calculated from center co-ordinate, z_c, as follows:

$$\begin{array}{l}{z}_{a1}=z_c-{\mathrm{e}}^{i\gamma }\cdot (0.5L+0.5W\cdot e^{i\beta })\\ z_{a2}=z_c-{\mathrm{e}}^{i\gamma }\cdot (0.5L-0.5W\cdot e^{i\beta })\\ z_{b1}=z_c-{\mathrm{e}}^{i\gamma }\cdot (-0.5L+0.5W\cdot e^{i\beta })\\ z_{b2}=z_c-{\mathrm{e}}^{i\gamma }\cdot (-0.5L-0.5W\cdot e^{i\beta })\end{array}(\mathrm{m})$$

(9b)

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), $\varphi$, and the stream function (Figure 15c), $\psi$, are mapped by the real part and the imaginary part of the complex potential respectively (see Equation (1)).

3.3. 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) which results in branch cuts in the potential function solution (Figure 15a) of Equation (9c), manifested as discontinuities in the pressure field when solved with MATLAB. The branch cuts do not appear in the stream function solution (Figure 15c) based on Equation (9d).

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 Figure 3b, Figure 4a, Figure 10a and Figure 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–d. However, the flow in Figure 15b is from bottom to top along the streamlines converging in the central flow domain, unimpeded by any adverse pressure gradients.

3.4. 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 $\mathsf{\Omega }(z)=\varphi +i\psi$ 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.

4. 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 Section 4.1 and Section 4.2, respectively. Section 4.3 and Section 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.

4.1. 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, $\mathsf{\Re }[\mathsf{\Omega }(z,t)]$. The complex potential is split in its real and imaginary parts, which correspond to respectively the potential function, $\varphi (z,t)$, and the stream function, $\psi (z,t)$ as shown in Equation (1). Subsequently, the potential function, $\varphi (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]:

$$\mathsf{\Delta }P(z,t)=-\frac{\varphi (z,t)\eta }{k}(\mathrm{Pa})$$

(10a)

where $\varphi (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² nor Darcy). The generic form of Equation (10a) 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, P_0, of the reservoir is known:

$$P(z,t)=P_0+\mathsf{\Delta }P(z,t)=P_0-\frac{\varphi (z,t)\eta }{k}(\mathrm{Pa})$$

(10b)

The porosity may appear as a scaling factor accounted for in input values for $\varphi (z,t)$ such as a certain flow strength [18].

4.2. Particle Tracking and Time-of-Flight Contours

The algorithm used to trace streamline and time-of-flight contours apply a Eulerian particle tracking method ($z_{n+1}\approx z_n+v(z_n)$), with the required velocity field inputs obtained from the velocity potential V(z) that can be derived from the complex potential $\mathsf{\Omega }(z,t)$ as follows:

$$V(z,t)=\frac{\overline{d\mathsf{\Omega }(z,t)}}{dz}=v_x+i{v}_y({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(11)

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 z₀ at time t₀ = 0 and calculating the initial velocity. By choosing an appropriate time-step, Δt, the position $z_1(t_1)$ of the tracer at time t₁ = t₀ + Δt is:

$$z_1(t_1)\approx z_0(t_0)+v(z_0(t_0))\mathsf{\Delta }t(\mathrm{m})$$

(12a)

The position $z_j(t_j)$ of the tracer at any other time t_j is:

$$z_j(t_j)\approx z_{j-1}(t_{j-1})+v(z_{j-1}(t_{j-1}))\mathsf{\Delta }t(\mathrm{m})$$

(12b)

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

4.3. 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⁻¹ (4.75 × 10⁻⁸ m·s⁻¹) 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⁻¹⁵ m²). 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⁻¹, 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). Figure 17b,c (left column) show the effect of flow accelerators, modeled by areal doublet, on both the streamlines and the TOFCs. The flow accelerators (which may mimic natural fractures) (Figure 17b,c) locally enhance the flow rate and some TOFCs reach further downstream, up to 37.6 m and 38.4 m in Figure 17b,c, respectively. The corresponding pressure plots are shown in Figure 17a–c (right column). Comparison of the pressure plots for each case shows that the presence of enhanced permeability flow channels increases the pressure gradient close to the flow elements. The inputs used to generate the Figure 17a–c are presented in Appendix B (

Table A2. Inputs for streamline simulation and pressure gradient calculation.

Inputs for Cases A through E
Flow channel permeability (k)	9.87 × 10−15 m2(10 mD)
Matrix Porosity	0.1
Viscosity (µ)	0.001 Pa.s (1 cP)
Simulation time for Euler’s particle tracking	6.31 × 108 s (20 years)
Time steps for particle tracking (Δt)	8640 s (0.1 days)

through

Table A4. Attributes of idealized synthetic flow channels (Case B).

Location	Width (m)	Length (m)	Angle (°)	Strength (m4·s−1)		Far-Field Velocity (m·s−1)
20	1	5	45	2.38	×10−6	4.75 × 10−9
20	1	5	315	2.38	×10−6

).

4.4. Effect of Spatial Complexity in the Location of the Flow Channels (Cases C–D)

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⁻¹ (or 4.7 × 10⁻⁸ m·s⁻¹), which maintains a maximum pressure of 2.9 × 10² kPa for a reservoir with a porosity of 10%. Figure 18b,c 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⁻⁸ and 9.4 × 10⁻⁶ (m⁴·s⁻¹). Figure 18b,c show that the pressure contour patterns become complex with addition of flow channels with varying strengths and orientations.

4.5. Field Example from Aravalli Supergroup—Rajasthan (India) (Case E)

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⁻⁹ m s⁻¹ 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. Figure 20b–d assume strengths of the flow channels with unit length scaled by one, two and five times 4.75 × 10⁻⁸ m⁴s⁻¹. Each flow channel is assumed to have a width of 0.01 m and unit height. The time-of-flight contours (Figure 20b) show that the fluid velocity increases near the regions with a high density of flow channels. The fluid which are not close to the flow channels travels through the matrix unchanged. For the cases with higher fracture strength such as in Figure 20c,d, the majority of the fluid will travel through the conduits and travels further downstream in the same run time.

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 high-pressure 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 (Figure 18b,c and Figure 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).

5. Discussion

5.1. 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 (Figure 1 and Figure 2) and the singularity doublet (Figure 7 and Figure 8). For the far-field flow the invariance of flows in any direction is apparent from:

$$\mathsf{\Omega }(z)=u_{\infty }{e}^{-i\gamma }z({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(13)

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 Section 2.2.3, Section 3.1 and Section 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 Figure 15b and Figure 16c, are invariant. Other placement of branch cuts (see later in Section 5.2) would render the result variant.

5.2. 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 multi-valued, 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 single-valued 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 (Figure 15a and Figure 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).

5.3. 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 closed-form 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 (Figure 17 and Figure 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 single-valued results (Figure 17, Figure 18, and Figure 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.

6. 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 (Figure 17 and Figure 20). Streamline and TOFC plots and pressure field solutions (Figure 17, Figure 18, and Figure 20) show that the presence of such flow channels have a significant impact on the flow in porous media. The algorithms and examples given in our study provide an efficient method to model the flow and pressure fields in fractured reservoirs.