## Methodology

Our modeling approach using CAM, side-steps the drawbacks of LGR and other gridded solutions (outlined in Section 2.2), by using a grid-less and meshless solution strategy, while still being able to model transient flow by time-stepping streamline integrals with very small time intervals to achieve smooth particle paths [74]. Single-phase flow in porous media can be modeled with analytical descriptions as has been advocated in many classical studies (e.g., Muskat [75,76]). The closed-form solution is commonly used to map curvilinear line-integrals to identify streamlines [77,78].

Our novel approach using complex analysis methods has been applied to fluid flooding and sweeping of hydrocarbon reservoirs with finite boundaries [6] and infinite boundaries [79,80]. Investigations of how hydraulic fractures drain the matrix in shale reservoirs are also possible using our method for streamline simulation and fluid time-of-flight tracing based on complex analysis methods.

Applications include the visualization of hydrocarbon migration toward hydraulically fractured wells [9], and fluid drainage near multi-fractured horizontal wells with fracture hits [8,74,81]. The key algorithms and key workflow steps used in our model method are explained below.

### Complex Potential and Generic Velocity Potential

The complex potential Ω(z) links the stream function (ψ) and potential function ф]):

The velocity components (v_{x} and v_{y}) can be obtained using [v_{x}=dф/dx=dψ/dy] and [v_{y}=dф/dy=-dψ/dx]. The full velocity field V(z) is given by the velocity potential:

The fluid flow is assumed irrotational, incompressible, immiscible and isothermal, such that the viscosity and density of the reservoir fluid remains constant. Capillary pressure, wettability and relative permeability effects are neglected. Details on complex analysis and flow applications can be found elsewhere [82,83,84,85,86,87].

### Single Point Source

The interval source and dipole elements used in our study are based on the superposition of point sources and point sinks. The complex potential centered of a source/sink flow centered at the origin with time-dependent strength m(t), is:

The flow in a porous medium due to a pressure gradient near a vertical injector wellbore in the planar (x,y) field is modeled by a point source (Figure 2a). Likewise, flow toward a vertical producer wellbore is represented by a point sink (Figure 2b). The corresponding velocity field is given by [79,80]:

where m(t) is the well strength (producers m_{(t)} > 0; injectors m_{(t)} < 0) in wells located at z_{c}.

**Figure 2.** (a) Point source, (b) point sink. Adapted from Weijermars and Van Harmelen [79].

Strength m_{(t)} for a well with productivity q(t) (m^{3}·s^{−1}), reservoir height h (m), and reservoir porosity n is:

Factor B is the volume factor defined by the ratio of any oil volumes under reservoir conditions divided by that oil volume in the stock tank at surface conditions.

**Figure 3.** Point source injection with associated streamlines (blue) and time-of-flight (TOF) contours (red) spaced for 1 year; total TOF is 10 years. Parameters: zc = 0 + 0 × 1i, q(t) = 1.84 × 10−7 m3/s, h = 1 m, n = 0.2, and Δt = 0.1 day. After Van Harmelen and Weijermars [9].

Figure 3 gives streamlines and time-of-flight contours for a homogenous horizontal reservoir due to fluid injection from a single vertical well with constant injection rate. Radial flow results in isochronous time-of-flight contours that are circles.

#### Interval Source

The hydraulic fractures in our model were modeled by interval-sources, obtained by superposing an infinite number of point sources (Figure 4). Assuming an interval source centred at the origin, with total length L and strength m(t), the complex potential is [7]:

**Figure 4.** Flow superposition principle diagram. (a1) Point source array, (a2) Point sink array, (b1) Linear interval source, (b2) Linear interval sink. Adapted from Weijermars and Van Harmelen [79].

Rotation of the interval-source by angle β and shifting its centre from the origin to coordinate z_{c} (Figure 5) yields:

**Figure 5.** Linear interval sink representing fracture element centered on zc, with end-points za and zb, total length L, and tilt angle β.

The complex potential for N interval-sources, each with their own angle β_{k}, centre z_{ck}, length Lk and strength m_{k}, reads [8]:

Differentiation of the above expression with respect to z, results in the velocity field expression for N interval sources with arbitrary orientations [8]:

Streamlines and time-of-flight contours for the area drained by a vertical fracture are obtained with our streamline tracing algorithm and visualized in Figure 6. The isochronous time-of-flight contours are confocal ellipses.

**Figure 6.** Fluid injected via hydraulic fractures of uniform conductivity with streamlines (blue) and time-of-flight contours (TOFCs, red). Total TOF is 10 years and TOFC spacing is 1 year. Interval source parameters: z_{c} = 0 + 0 × 1i, q(t) = 1.84 × 10^{−7} m3/s, h = 1 m, n = 0.2, L = 5 m, β = 0° and Δt = 0.1 day. After Van Harmelen and Weijermars [9].

### Single Point Dipole

The micro-cracks in our study were modeled by line dipoles. The point dipole, also known as point doublet, singularity doublet and singularity dipole, is a common occurrence in amongst other electromagnetism [88], aerodynamics [89], and fluid mechanics [90]. Deriving the point dipole velocity field entails combining the velocity fields of a point source and point sink (Figure 7) in a limiting process. Placing a point sink and point source of equal but opposing strengths on the real axis at a distance of 2ε from each other, meaning their locations are ε and −ε respectively, yields the velocity field:

**Figure 7.** (a) Point sources and (b) point sink, superposed to form (c) singularity dipole (adapted from Weijermars and Van Harmelen [79]).

The point dipole is then obtained by letting 2ε approach zero while the product q(t)∙2ε remains constant (meaning q(t) increases inversely proportional to distance 2ε). By defining this constant as υ(t) = q(t)∙2ε, we find that the velocity field for the point dipole, located at the origin, is:

Rotating the point dipole of expression (11) counter-clockwise by angle β and placing it at location zk is possible by employing the conformal mapping f(z) = e^{−βi}(z − z_{k}) and evaluating V(f(z))·f’(z). This results in the velocity field formula:

Note that υ is the strength of the point dipole and, as a result of the limiting process, its unit is (m^{4}·s^{−1}). The complex potential of a point dipole is hence given by:

### Line-Dipole

The analytical point dipole in expression (12) has its source side pointing parallel to the real axis, to the right, if β = 0 (Figure 8). We derive a line-dipole, starting with β = 0 and its center at real axis coordinate xc. The resulting velocity field is:

**Figure 8.** (a) Array of singularity dipoles, and (b) line dipole (adapted from Weijermars and Van Harmelen [79]).

Next combine j of the point dipoles, with a uniform strength υ(t):

To transform this expression into a Riemann integral, the interval [−0.5L, 0.5L] is partitioned into j intervals, by defining the j + 1 points:

Therefore, each point dipole is located at:

In order to turn expression (15) into a Riemann integral, the spacing between two point dipoles (Δx_{k}) needs to be defined:

Multiplying expression (15) with L/L and splitting the sum then gives:

By letting j increase without bound, the last term in expression (19) goes to zero and a Riemann integral is obtained because Δxk goes to zero. The resulting integral and velocity field is:

Rotating this expression and shifting the center is again achieved with the conformal mapping f(z) = e−β_{i}(z − z_{c}) and by evaluating V(f(z))·f′(z). The velocity field for the line-dipole is therefore:

Expression (21) is similar to the one in Weijermars and Van Harmelen [79], though more explicitly formulated. The complex potential of a line-dipole is hence given by:

### History Matching and Flux Allocation

A general expression used for history matching of the prototype well rate is given by Duong [91] (using field units in bbls/day):

Cumulative oil production is:

The next step is to allocate total well production rates to the individual fracs in our drainage model:

C = 5.61458333 converts field input units of [Math Processing Error] in stb/day to field output units of [Math Processing Error] in ft^{3}/day, and WOR = Swater/Soil, the respective saturations of water and oil.

In our study, the total type curve output, q_{well}(t), was simply prorated to the stages studied (13 out of 33) and divided proportional to the surface area of each frac stage (as determined by the product of frac heights, h_{k}, and lengths L_{k}) relative to the total surface area of the 13 fracs, and correcting for the water to oil ratio, WOR:

The prorated flux per frac over time was used in the drainage visualization of our study.

### Oil in Place and Recovery Factor

The produced oil can be compared to the original oil in place:

With field units inputs used as follows: 7758 = barrels (bbls) of oil per (acre·ft); A = area in acres (TBD in acres); h = height in feet; n = porosity; Sw = water saturation; B = volume factor.

The oil recovery factor is:

with F = Recovery factor, and Np(t) the cumulative production (stb) till a certain time t.