Geofluids
Volume 2017 (2017), Article ID 5609571, 18 pages
https://doi.org/10.1155/2017/5609571

Research Article

Rules for Flight Paths and Time of Flight for Flows in Porous Media with Heterogeneous Permeability and Porosity

Lihua Zuo and Ruud Weijermars

Harold Vance Department of Petroleum Engineering, Texas A&M University, 3116 TAMU, College Station, TX 77843-3116, USA

Correspondence should be addressed to Ruud Weijermars; r.weijermars@tamu.edu

Received 24 April 2017; Revised 18 November 2017; Accepted 3 December 2017; Published 31 December 2017

Academic Editor: Marco Petitta

Copyright © 2017 Lihua Zuo and Ruud Weijermars. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Porous media like hydrocarbon reservoirs may be composed of a wide variety of rocks with different porosity and permeability. Our study shows in algorithms and in synthetic numerical simulations that the flow pattern of any particular porous medium, assuming constant fluid properties and standardized boundary and initial conditions, is not affected by any spatial porosity changes but will vary only according to spatial permeability changes. In contrast, the time of flight along the streamline will be affected by both the permeability and porosity, albeit in opposite directions. A theoretical framework is presented with evidence from flow visualizations. A series of strategically chosen streamline simulations, including systematic spatial variations of porosity and permeability, visualizes the respective effects on the flight path and time of flight. Two practical rules are formulated. Rule  1 states that an increase in permeability decreases the time of flight, whereas an increase in porosity increases the time of flight. Rule  2 states that the permeability uniquely controls the flight path of fluid flow in porous media; local porosity variations do not affect the streamline path. The two rules are essential for understanding fluid transport mechanisms, and their rigorous validation therefore is merited.

1. Introduction

Flow analysis in porous media is at the macroscopic scale governed by Darcy’s law. The local flux and velocity are controlled by reservoir parameters such as permeability and porosity and by fluid properties such as density and viscosity. However, permeability and porosity are not directly connected, neither by their formal definition nor by dimensional units. In fact, porosity [fraction of void space] is a static material property of a porous medium. Permeability [m²] is an independent dynamic scaling parameter in Darcy’s law that relates the fluid transmission flux [m⁻³ s¹] through a unit area [m²] by the ratio of the applied pressure gradient [Pa m⁻¹] and the transmitted fluid’s dynamic viscosity [Pa s].

Although physically independent, there often exists an empirical relationship between the permeability and porosity for a given porous medium [1–4]. The well-known Carman-Kozeny relationship links the porosity and permeability in nondimensional form using a characteristic flow tube dimension for scaling [5–7]. The Carman-Kozeny relationship is a theoretical model which has limited practical value for reservoirs comprised of rock types with exceptional high porosity and relatively low permeability [8, 9]. For example, Figure 1(a) shows a microscopic example of a diatomite carbonate, where relatively high porosity (~60%) combines with very low permeability (~1 mD). In practice, reservoirs properties are often dominated by certain rock types with petrophysical properties that correlate permeability and porosity in certain domains (Figure 1(b)).

(a)

(b)

Figure 1: (a) SEM image of porosity structure for diatomite (from Erdőbénye, Hungary, see [12]). (b) Domains of porosity and permeability correlation for unconsolidated sands and various rock types based on Nelson [2], revised, and adding point 1(a) for diatomite sample (the permeability is in log scale).

Darcy’s law features permeability, but porosity does not directly appear. However, the equation of motion includes both permeability and porosity, which are important parameters for predicting, respectively, the flow paths (FP) and time of flight (TOF) of fluids transporting in porous media [10, 11]. The fluid properties may also complicate any flow analysis in porous media when temperature and pressure effects change the viscosity, density, and phase behavior. However, if we assume constant viscosity and single phase, incompressible flow, then the fluid properties may only cause pressure variations, but these do not affect any flow path nor rate of flow when wells are pumped at constant rates. In a reservoir space with single phase, isoviscous and incompressible flow, not the fluid properties, but only the permeability and porosity control the flow paths and time of flight between any array of injector wells and producer wells. We will elaborate the details of how permeability and porosity relate to FP and TOF in the remainder of this paper.

Although for any particular rock type a higher permeability generally correlates with higher porosity (Figure 1(b)), these parameters each independently affect TOF in opposite directions. The time of flight changes with the porosity: flow speeds up when the porosity decreases (thus shortens TOF) and flow slows down when the porosity increases (thus lengthens TOF). For permeability the opposite occurs: flow speeds up when the permeability increases and slows down when the permeability decreases. The corresponding effects on TOF are indicated in Figure 1(b) along the permeability and porosity axes.

There is considerable risk for drawing overly quick conclusions in practical situations. Consider a simple doublet composed of one injector and one producer (Figure 2(a)) with water breakthrough much faster than anticipated (actually twice as fast, Figure 2(b)): there are several possible explanations. One conclusion is that the average reservoir permeability is higher (actually factor 2) than initially assumed. However, an alternative, equally valid interpretation could be that the permeability assumption was correct, but the porosity, which also affects the time of flight and flood arrival time, is only 1/2 of the initially assumed value. A third explanation would be that the permeability was 50% underestimated and porosity 25% overestimated.

(a) A priori model

(b) Posterior observation

Figure 2: Example  1: Illustration of possible misinterpretation of faster water breakthrough. Streamlines are tracked in blue and red curves are the time of flight contours spaced for 0.25 years. Water breakthrough takes 1 year in a priori model (a) but the posterior observation shows breakthrough occurs already after half a year (b).

In another example, consider a localized zone of faster flow, such as a longitudinal hydraulic fracture connecting two vertical wells in Figure 3. The streamlines are concentrated in the high conductivity fracture, which could be easily interpreted as a consequence of increases in both the porosity and permeability (due to the fracturing process). However, careful assessment demonstrates that streamlines bungle into an hydraulic fracture only because of the local increase in permeability. The porosity change does not affect the streamline path. Any porosity increase would merely slow the flow rate in the hydraulic fracture but is not responsible for any divergence or convergence of the streamlines.

Figure 3: Example  2: Streamline trajectories passing through a hydraulic fracture zone. Flow plane is perpendicular to wells and fracture surface.

Our study intends to deliver a practical set of rules for the relationships highlighted in Figures 2 and 3, with explicit theoretical proof as well as numerical experimental validation. In our opinion, a combination of theoretical derivation with equations and visualization for specific cases with numerical coding is essential. We illustrate these fundamental aspects of fluid flow in porous media with heterogeneous permeability and porosity (for both isotropic and anisotropic cases), using systematic flow visualizations based on the integration of three software packages, calibrated with an independent analytical flow description. In addition to the general rationale given above, the Discussion provides examples of laboratory studies of permeability anisotropy where permeability is the main focus and porosity is not adequately accounted for.

2. Research Methodology

2.1. Reservoir and Well Design

The basic flow visualization model uses 5 injectors and 5 producers in a direct line drive arrangement (Figure 4) geometrically similar to that used in our earlier study benchmarking an analytical model [13]. The well flow rates are fixed over time so that the flow is in the steady state, in which case the streamlines are coincident with the pathlines of the flows. For the simplicity and without loss of generality, we assume the following properties for the reservoir. Lateral flow space is relatively large and we assume a thin reservoir with no vertical velocity gradients. In essence we have a 2D flow space but can account for volumetric flux by assuming unit thickness of the reservoir space. Lateral flow space is taken large enough to avoid any boundary effects on the streamline patterns. Fluid flow occurs in piston-like displacement of the oil-water contact so that we study a single-phase flow problem (Figure 4). This is required in order to be able to focus exclusively on the effects on the FP and TOF due to local variations in the permeability, porosity, and anisotropy of the porous medium and of any fluid viscosity changes.

Figure 4: Well constellation used in this study assumes one horizontal injection well with water injection rates controlled by 5 inlet control valves (ICVs). Flooding occurs by direct line drive between the water ICVs and 5 vertical producer wells where pressure can be monitored with installed bottomhole assemblies (BHAs; after [13]).

Previous theoretical work and field experiments have demonstrated that the relationship between injected fluid and production is not always simple [14, 15]. The basic merit of streamline simulations is that well productivity can be explained as the flux of fluid carried through the streamtubes outlined by discrete bundles of streamlines into the well (Morel-Seytoux 1965 [16–21]). An excellent review of streamtube reservoir models has been given by Thiele [22] and Thiele et al. [23]. Parsons [24] was the first to consider permeability anisotropy, which substantially alters recovery rates as compared to isotropic reservoirs. Directional peaks in the interconnectivity coefficients highlight high permeability channels in the reservoir and vice versa. An advantage of flow-based allocation [25] over allocation factors based on fixed streamtubes are so-called time of flight models [26], which account in the well capacitance sum for transient flow effects such as well shut-ins, infill drilling, changes in injection rate, and/or field pressure decline.

The porosity of porous media could be divided into connected porosity and unconnected porosity. Connected porosity could be defined as the ratio of the volume of fluid that can flow into the rock, while fluids cannot access unconnected pores. In this study, we assume the porosity is connected porosity. In our present study, we avoid any transient effects by imposing stationary, creeping flow and assume incompressible fluids. Under this assumption, the flow rates of injectors and producers will be kept constant. FP are fixated by the particle paths for transient flows and by streamlines for steady flows. Our work provides such proof and validation for the fact that the TOF is controlled by the integrated velocities of fluid particles as they move along the flight path.

2.2. Scale of Study

In a porous medium, the material which is homogeneous at the scale of 1 millimeter might present heterogeneities at the scale of 100 meters. Unconventional reservoirs require dealing with multiscale problems due to the very low permeability in shale gas and oil reservoirs. On the microscopic scale, the fluid may perform as non-Darcy flow where classical Darcy’s law will no longer stand. To better model non-Darcy flow, different tools and flow equations have been proposed to simulate the atomic level of the diffusion phenomena; for example, the fractional diffusion equations were proposed to replace the classical diffusion equation [27–32]. The fractional decline curve model was introduced for well rate analysis and production forecast [33, 34].

While new tools might be useful for some unconventional reservoir plays, currently, for most popular models and simulators in petroleum industry, the conventional flow equations such as the stationary Stokes flow including Darcy’s law are still widely used. These equations are generally built on continuum mechanics, which is based on the assumption that bulk properties such as permeability, density, and porosity, all, vary continuously over the continuum on the macroscopic scale (usually between millimeter to kilometer). In continuum mechanics, at each point of the continuum, each specific bulk property is represented by the average over a representative elementary volume [35]. In the present study, we work on the macroscopic scale and assume that Darcy’s law stands.

2.3. Theoretical Derivation

We will show below the theoretical proof of our contention that porosity does not influence the flight path outlined by the streamlines, whereas porosity and permeability both influence time of flight. Such insight follows from the two rules postulated in our Abstract: Rule  1 states that an increase in permeability decreases the time of flight, whereas an increase in porosity increases the time of flight. Rule  2 states that the permeability uniquely controls the flight path of fluid flow in porous media; local porosity variations do not affect the streamline path.

We will first deliver explicit theoretical proof for the above rules. The physical parameters that control the streamline trajectories and the time of flight are solely due to reservoir properties and initial conditions. First, consider the following Stokes equation for stationary, creeping flow of an incompressible fluid [35]:

where is the viscosity, is the velocity in the direction, is the gravity component in the direction, and is the pressure. For a porous medium, the viscous resisting force is not only affected by the viscosity and the pressure gradient but also affected by properties of the porous medium captured in porosity, , and the permeability tensor, :

To express which parameters control the velocity in the direction, multiply (1b) by on both sides and divide by :

From (1c), we can see that for a certain pressure gradient the fluid velocity is linearly related to permeability and inversely proportional to both porosity and fluid viscosity.

If next we multiply with both sides of (1c), we obtain the volumetric flux in the direction:

In isotropic porous media the off-diagonal elements in the permeability tensor are zero and the diagonal elements are identical, so that we arrive at common Darcy’s law:

Correspondingly, the velocity vector due to a given pressure gradient is given by

which is similar to (1c). Equation (2b) shows that the velocity is linearly proportional to the permeability and inversely proportional to the porosity, which is in accordance with Rule  1.

To show Rule  1 more directly, the relation between permeability/porosity and TOF is derived next. It is known that the time of flight, , is defined as [11]

where is the spatial distance along streamline. Combining with (2b), we can see that for a certain pressure gradient and a given fluid viscosity the time of flight is related to both permeability and porosity:

Again we can reason spatial changes in either the porosity or viscosity will cause only a proportional change in , due to which expression (3b) can be simplified into (ignoring gravity effects)

with constant accounting for a given combination of pressure gradient, porosity, and viscosity. Equation (3c) shows directly that TOF is linearly related to the porosity (faster arrival due to faster flow when decreases). Time of flight is shorter for higher permeability, with an inverse linear relationship. Consequently, porosity and permeability have opposite effects on TOF (Rule  1).

An alternative expression relates the velocity field to the stream function; in the 2D case, we have the following corresponding relation:

For a given streamline pattern to remain constant throughout the 2D flow space, the absolute velocities may change but the ratio of two velocity vector components may not change.

Combining expressions (2b) and (3d) gives the following relation between the velocity vector components and any spatial variations in properties of the reservoir:

We consider a 2D problem in a horizontal plane and assume gravity components do not spatially change on the scale of the flow and thus can be neglected (i.e., ). Note that is isotropic in our current study. Fluid viscosity remains constant and incompressible flow implies special density changes cannot occur: = constant. Combining expressions (4a) and (4b) with the given assumptions results in the following relation between the spatial streamline tangents and properties of reservoir:

For a given pressure field any spatial changes in the porosity will cause only a proportional change in (see the Appendix for more explanations of this); the ratio will remain constant, , due to which (4c) may be simplified into

Equation (4d) shows that the velocity field and implied particle paths (streamlines for steady flows) will be affected solely by spatial changes in permeability; there is no relation with porosity. This supports Rule  2.

After presenting and demonstrating the theoretical development in the previous paragraphs, next we will use numerical experiments to validate the inferred relationships between streamline trajectory patterns and permeability [see (4d)] and the controlling factors of time of flight (porosity and permeability [see (3c)]). The advantage of our systematic reservoir simulations is that the streamline patterns and time of flight contours could be visualized in a systematic fashion. Without such visualization of the streamline trajectories under controlled input parameters, one cannot easily verify whether or not a streamline path is affected by local changes in the porosity. A similar benefit occurs when verifying the impact on the time of flight due to spatial changes in porosity and permeability. Additionally, we also show that for a particular constant pressure profile in the reservoir the velocity field is independent of fluid viscosity, but the pressure field will change by variations in porosity, permeability, and/or viscosity (see the Appendix).

2.4. Numerical Model Design

Our flow visualization experiments use a numerical streamline tracing algorithm, which allows comparison of time of flight as well as streamline patterns for the advancing flood front as it becomes affected by variations in reservoir properties. The base data used for our experiments are listed in Table 1.

Table 1: Numerical model parameters.

Parameters
Dimensions
Number of grids
Porosity
Permeability
Injectors
Producers

Value
250 × 200 × 1
250 × 200 × 1
20
100
Water, 1000
Oil, 1000

Unit
ft
ft

mD
bbls/d
bbls/d

The numerical streamline algorithm calculates the TOF information by using the methodology described in Datta-Gupta and King [11] and Zuo et al. [33, 34], which is an extension of Pollock’s [41] solution. Based on the data provided in Table 1, we first designed the model in ECLIPSE to obtain reservoir results such as pressure and flow rates on all six faces for each cell. Next we run the streamline algorithm to obtain the streamline tracing results, which are then imported in Petrel to visualize the streamlines. The numerical streamline tracing method is used below to systematically investigate and visualize the impact of porosity and permeability distributions on TOF and FP in porous media with heterogeneous permeability and porosity. The basic well design uses 5 injectors and 5 producers in a direct line drive arrangement geometrically similar to that used in our earlier study benchmarking and analytical model [13, 42]. We restrict our analysis to 2D visualization of water advancement in a horizontal reservoir layer as portrayed in Figure 4.

2.5. Numerical Model Verification

For simple reservoir models especially if it is homogeneous porous medium and simple injectors and producers placement as shown in Figure 4, an alternative, analytical streamline method can be applied to compute FP and TOF [13, 42–45]. The analytical streamline method is based on complex potentials and calculates the flow velocity using the strength/locations of the point source (injectors) and sinks (producers) by the following formula [13]:

where is the strength of each point source/sink and is the location of each source/sink. Using (5) combined with a first-order Eulerian scheme , the streamline trajectories and time of flight contours can be calculated.

The analytical streamline method and the numerical streamline algorithm have been cross-validated in our previous studies [13, 42]. The numerical streamline package used in the present study integrates sophisticated functionality of advanced reservoir simulators, such as a capacity to account for heterogeneities and other complex geological properties of the reservoir. Such a simulator is well suited for testing the impact on FP and TOF of variations of the permeability and porosity. For completeness of the present study, it includes comparisons between the two methods. Figures 5(a) and 5(b) are the analytical streamline results and numerical streamline results for a homogeneous reservoir, respectively. Figures 5(c) and 5(d) are the corresponding results for a reservoir with heterogeneous permeability distribution, namely, a higher permeability distribution (1000 mD versus 100 mD) in the regions between I4/I5 and P4/P5. Both the streamline trajectories (FP) and TOF are matched very well (see [13] for more details about the matches for the streamline trajectories and time of flight contours).

(a)

(b)

(c)

(d)

Figure 5: Verification of the numerical streamline method with the analytical streamline method (after [13]). (a) The streamline result (blue and yellow lines) and time of flight contours (red curves) of the analytical method for the base model. (b) The streamline result (red lines) and time of flight contours (rainbow color) of the numerical method for base model. (c) The streamline result and time of flight contours of the analytical method for the high permeability zone model. (d) The streamline result and time of flight contours of the numerical method for the high permeability zone model.

2.6. Model Design for Numerical Experimental Validation

The aim of our study is to investigate the effect of permeability and porosity on FP and TOF by some systematically designed numerical experiments. From the synthetic simulations, we can check our assertion that an increase in permeability decreases the time of flight, whereas an increase in porosity increases the time of flight (Rule  1). The permeability uniquely controls the flight path of fluid flow in porous media; local porosity variations do not affect the streamline path (Rule  2).

The first series of experiments is comprised of heterogeneous reservoirs including isotropic domains with both permeability and porosity contrasts.

 Case A is a reservoir with homogeneous permeability and homogenous porosity distribution. Case B is a reservoir with homogeneous permeability and heterogeneous porosity distribution. Case C is a reservoir with minor domains of lower permeability in a higher permeability matrix: matrix porosity is held constant but porosity of lower permeability domains is varied to show that in such low permeability domains TOF can still equal that of the higher permeability domains. Case D is a reservoir with minor domains of higher permeability in a lower permeability matrix: matrix porosity is held constant but porosity of the higher permeability domains is varied to show impact in TOF in such low permeability domains can still equal that of the lower permeability domains.

The second series of experiments is comprised of streamline tracing cases which include directional anisotropy of the permeability, while keeping porosity constant and isotropic for all cases.

 Case E includes permeability anisotropy. Case F includes slanted permeability zone. Case G includes multilayered heterogeneous permeability zones. Case H includes multilayered heterogeneous porosity zones.

2.7. Homogeneous and Heterogeneous Isotropic Reservoirs

This section presents results of flow visualizations including domains with both permeability and porosity contrast. First a homogeneous base case is presented, with reservoir and design parameters specified in Table 1.

2.7.1. Case A: Homogeneous Permeability and Homogenous Porosity

The flow visualization of Figure 6(a) is for a base case reservoir with permeability and porosity both uniform throughout the flow space. Water injection occurs at the lower array of 5 injectors and production occurs at 5 producer wells horizontally aligned near the upper half of the image space. Streamlines (red curves) connect one or more pairs of the injector and producer wells. The advance of the frontal edge of the waterflood is outlined by the yellow TOF contour. A further subdivision of the time of flight contours for locations passed by the waterflood at earlier stages is given by the rainbow colors, for which a TOF scale is given for each simulated case (Figures 6(a)–6(c)).

(a)

(b)

(c)

Figure 6: Top view of flow in horizontal reservoir with finite thickness but free flow between top and bottom boundaries (resp., above and below the image plane). Horizontal spacing between wells is 6 ft and vertical separation of injector and producer well arrays is 30 ft. Streamlines (red curves) and flood flight times (rainbow colors) are based on ECLIPSE pressure data augmented with proprietary velocity calculations based on the flux balance of adjacent cell nodes. Scale bar gives TOF in days and is valid for all three cases (a–c). Water front advance is outlined by the yellow leading edge as the waterflood moves away from the 5 injectors (lower row). (a) Base case permeability 100 mD and porosity 10%, TOF = 0.055 days. (b) Increasing porosity to 20% results in linear increase of TOF = 0.11 days. (c) Further increase of porosity to 40% results in linear increase of TOF = 0.22 days.

Figure 6(b) keeps the permeability unchanged but doubles the porosity from 10% to 20%, due to which the time of flight is slowed by half and TOF doubles (Figure 6(b)). However, the streamline pattern of Figure 6(b) is indistinguishable from that shown in Figure 3(a). Further doubling of the porosity from 20 to 40% confirms the linear relationship between TOF and porosity (Figure 6(c)). Note that streamline paths for Figures 6(a)–6(c) remain identical, in spite of porosity changes.

2.7.2. Case B: Homogeneous Permeability and Heterogeneous Porosity Distribution

We next keep the permeability constant throughout the reservoir but the porosity is assumed to be heterogeneous, made up of two distinct domains that have a porosity, respectively, lower and higher than the ambient reservoir. The domain of lower porosity results in a local shortening of TOF (Figure 7, yellow box) due to speeding of the flow in the lower porosity domain. The domain of higher porosity results in a local increase of TOF (Figure 7, blue box) due to slowing of the flow in the higher porosity domain. In neither case is the streamline pattern affected; the path of the streamlines is entirely determined by the permeability distribution alone. This is in agreement with the Law of Streamline Refraction, which says streamlines are only affected when fluid crosses boundaries with permeability contrast [35]. Our flow simulations visualize conclusively that porosity variations do not affect the flight path outlined by the streamlines. However, the time of flight will vary proportionally with porosity changes (Figure 7).

Figure 7: Matrix and domains of heterogeneous porosity both have base case permeability of 100 mD. Matrix has 40% porosity, except for domains of the yellow box (5% porosity) and blue box (80% porosity).

2.7.3. Case C: Minor Domain of Lower Permeability

The porosity is held constant throughout the reservoir except for one low permeability domain (rectangular region) between the injector and producer wells where the porosity is systematically but uniformly varied to show the additional impact on TOF in such low permeability domains. For example, a high permeability zone will diverge streamlines such that more streamlines cross the higher permeability domain (Figure 8(a)). If the porosity remains the same throughout the porous medium (20% as in our base case), the slower fluid throughput in the low permeability zone results in longer TOF so that fluid moves across the zone with a decreased flux. Streamlines are diverging in the heterogeneous zone due to slowdown of fluid flow (Figure 8(a)). When the porosity in the low permeability domain is adjusted from 20% to 5%, the decrease of TOF in the low permeability domain is undone (Figure 8(b)). Decrease in the time of flight due to a lower permeability zone is countered when such zones have a relatively low porosity, which shortens the TOF by speeding up the fluid flow. The cause of the flow speed increase is that with unchanged permeability but lower porosity the same fluid flux needs to cross a smaller space. Less space decreases the TOF and the porosity drop can fully counter the effect of a lower permeability (Figure 8(b)). The streamline patterns of Figures 8(a) and 8(b) are the same; slight differences in appearance are due to different particle bundles being tracked in each flow simulation. Figure 8(c) illustrates that TOF can be increased (longer TOF) still further (as compared to the base case of Figure 8(a)) when the porosity in the heterogeneous region is adjusted from 20% to 80%. The larger porosity leads to a further increase in the TOF in spite of the lower permeability in the central heterogeneity.

(a)

(b)

(c)

Figure 8: Streamlines (red curves) and flood flight times (rainbow colors). (a) Domains of lower permeability result in local increase in TOF. Matrix permeability is 100 mD; rectangular central zone between injectors and producers has permeability of 20 mD. Porosity is equal in both domains (20%). (b) Increase of TOF in low permeability domain is undone when porosity in low permeability domain is adjusted from 20% to 5%. (c) Reversely, TOF can be increased still further than in (a) when porosity is adjusted from 20% to 80%.

2.7.4. Case D: Minor Domain of Higher Permeability

Figures 8(a)–8(c) illustrated the effect of a relatively low permeability heterogeneous zone set in a higher permeability matrix. We now consider the reverse case of a relatively high permeability heterogeneous domain surrounded by a less permeable matrix. Streamlines are converging in the heterogeneous zone due to speeding up of fluid flow (Figure 9(a)). The TOF is reduced in the high permeability zone when the effective porosities of the heterogeneous domain and ambient reservoir are equal. Next, matrix porosity is held constant but porosity of the high permeability domain is varied to show impact on TOF of porosity changes in such higher permeability domains. When the porosity in the low permeability domain is adjusted from 20% to 80%, the increase of TOF in the high permeability domain is undone (Figure 9(b)). Increase of the time of flight due to a high permeability zone is countered when such zones have a relatively high porosity, which increases the TOF by slowing down fluid flow. The reason of flow speed decrease is that with unchanged permeability but higher porosity the same fluid flux needs to cross a larger space. More space slows down the TOF and the porosity change fully counters the effect of the larger permeability (Figure 9(b)). The streamline patterns of Figures 9(a) and 9(b) are the same; slight differences in appearance are due to different particle bundles being tracked in each flow simulation. Figure 9(c) illustrates that TOF can be fastened still further than in case Figure 9(a) when porosity is adjusted from 20% to 5%.

(a)

(b)

(c)

Figure 9: Streamlines (red curves) and flood flight times (rainbow colors). (a) Domains of higher permeability result in local decrease in TOF. Matrix permeability is 100 mD; rectangular central zone between injectors and producers has permeability of 1000 mD. Porosity is equal in both domains (20%). (b) Decrease of TOF in high permeability domain is undone when porosity in low permeability domain is adjusted from 20% to 80%. (c) Reversely, TOF can be decreased still further as compared to (a) when porosity is adjusted from 20% to 5%.

2.8. Directional Permeability Anisotropy

This section considers results for reservoirs with directional permeability anisotropy while keeping porosity constant throughout the reservoir.

2.8.1. Case E: Permeability Anisotropy

It is well known from the early work of Muskat [46, 47] that streamlines will be affected by permeability anisotropy. Figures 10(a) and 10(b) visualize two cases of permeability anisotropy. Permeability anisotropy in the direction normal to the general flow direction of the direct line drive will slow down the water flood due to outward migration of the streamlines as compared to the homogenous uniform permeability case of Figure 8. In contrast, permeability anisotropy with the higher permeability direction aligned with the general flow direction of the direct line drive will fasten the water flood advance with streamlines running closer and more parallel to each other compared to the base cases of Figure 8. Water flood design strategies should take such effects into account when zones of predominant permeability anisotropy occur in the reservoir. More results and permeability anisotropy measurement results can be found in Table 2 in the Discussion.

Table 2: Permeability anisotropy studies.

Formation location
Horizontal permeability (nD)
Vertical permeability (nD)
Ratio (horizontal/vertical)
Measuring method
Reference

Barnett
96.3
2.3
45
Pulse decay method
Bhandari et al., 2015 [36]

Eagle Ford
50–210
50–80
4
Pulse decay method
Mokhtari and Tutuncu, 2015 [37]

Wilcox
11800
500
23
Transient pulse method
Kwon et al., 2004 [38]

Longmaxi
199.3–201.3
8
25
Pulse decay method
Pan et al., 2015 [39]

Krechba
12–58501.4
0.6–264
10–100
(one is 50000)
Transient pulse method
Armitage et al., 2011 [40]

(a)

(b)

Figure 10: Impact of permeability anisotropy on FP and TOF. Streamlines (red curves) and flood flight times (rainbow colors). (a) Permeability in the horizontal -direction equals 100 mD while permeability in the vertical -direction equals 20 mD. (b) Permeability in the -direction equals 100 mD while permeability in the -direction equals 500 mD. Porosity is uniform and equal in both simulations (20%).

2.8.2. Case F: Slanted Inhomogeneous Permeability Zone

Local permeability anisotropy may occur in many situations. In clastic sediments, it may occur along fractures and fault zones, which may either increase permeability due to dilational effects, or reduce permeability, for example, due to clay smear. For example, carbonate-bearing fault zones can be interpreted as damage zones with the permeability increasing 2 orders of magnitude relative to that of the undamaged protolith [48]. In many other cases, the nonuniform distribution of deviatoric stress could also lead to local permeability anisotropy. The impact on the flow of a dilational fracture with an increased permeability is illustrated in Figure 11(a). The slanted high permeability zone refracts the streamlines. The TOF is shortened by the faster flow across the fracture zone. The rightmost producer sees first water breakthrough. Figure 11(b) illustrates the impact of a low permeability zone, causing streamlines to refract such that these align with the trend of the fault. The angular change in the streamline orientation is given by the Law of Streamline Refraction [49].

(a)

(b)

Figure 11: Streamlines (red curves) and flood flight times (rainbow colors). (a) Domains of lower permeability result in local increase in TOF. Matrix permeability is 100 mD; upright to bottom left slanted zone has permeability of 10 mD. Porosity is equal in both domains (20%). (b) Matrix permeability is 100 mD; upright to bottom left slanted zone has increased permeability of 1000 mD. Porosity is constant and equal in all domains (20%).

2.8.3. Case G: Multilayered Heterogeneous Permeability Zones with High Conductivity Gaps

Our final set of simulations is shown in Figures 12(a) and 12(b). In this example, we use the same base case as defined in Table 1, but we add four extra low permeability layers such (e.g., clay) as what is denoted by the gray color in Figures 12(a) and 12(b). We assume a horizontal higher permeability sheet (e.g., sand) with vertical flow barriers and baffles made up of low permeability clay walls that prevent fluid flow, except for certain gaps in the clay/shale wall, assuming dilation fractures filled with sand. The difference between Figures 12(a) and 12(b) is that we only have one gap in Figure 12(b) compared to two gaps in Figure 12(a). Alternatively, the model in Figure 12(b) may be interpreted to represent a pore scale effect, where the shortest streamline connections provide the fastest flight path. The corresponding streamline results are shown in Figures 12(c) and 12(d), respectively, from which we can see that more tortuous flight paths inherently have a longer time of flight. The time of flight for the left side of model of Figure 12(d) will be almost the same as that in Figure 12(c), but the right part will have longer flight times due to the missing gaps in the top clay wall. Essentially, the flux contribution from the three injectors on the right in Compartment 1 of Figure 12(d) cannot arrive at that gap in time. We have also prepared two animations of the advective flood advancement of Figures 12(c) and 12(d), which visualize the dynamic flood migration [provided as auxiliary material with this article].

(a)

(b)

(c)

(d)

(e)

(f)

Figure 12: Two central clay or shale walls obstruct flow, except for sand filled fractures that leave a high permeability conduit for the flood advance across the two central low permeability walls. Low permeability walls also occur above and below the well arrays. (a) There are two gaps on the second clay wall from the top. (b) There is only one gap on the second clay wall from the top. (c) The streamlines (red curves) and flood flight times (rainbow colors for TOF scaled in unit of days) for case (a), where sand permeability is 100 mD and clay is 100 nD; porosity everywhere equals 20%. (d) The streamlines and flood flight times for case (b). (e) Pressure map for case (a) (rainbow colors scaled for pressure in psi units). (f) Pressure map for case (b). Flow is in steady state, due to which pressures remain steady.

The corresponding pressure maps for the models of Figures 12(a) and 12(b) are given in Figures 12(e) and 12(f). The most intriguing aspect of these pressure maps is the large pressure difference between the three sand compartments or say the fast pressure drop across each clay wall. For example, in the model of Figure 12(e) the sand compartment near the injectors (Compartment 1) has a relatively uniform pressure of around 10,000 psi, while in the second sand section (Compartment 2) the pressure drops to about 5,500 psi, and the compartment near the producers (Compartment 3) has a relatively uniform pressure of around 500 psi. Without the presence of any compartments the pressure would drop approximately linearly with distance between the injectors and producers. Our new finding is that, with narrow flow conduits between high permeability compartments (Figures 12(a)–12(f)), the pressure gradient becomes concentrated in the narrow connections between the compartments. The implication for field development is that overpressure may occur in the lower compartment with uniform, high pressure, which could pose a drilling hazard when infill wells are planned.

The compartments may themselves have more or less uniform pressure distributions. Figure 12(f) illustrates that when flow connections of the central sand section (Compartment 2) to the adjacent high pressure compartment upstream (Compartment 1) are less restricted than to the adjacent low pressure compartment downstream (Compartment 3), the pressure in the central sand section may be partially closer to that of the adjacent high pressure section. The pressure in the central sand (Compartment 2) of Figure 12(f) ranges between 8,000 and 6,000 psi, whereas the pressure of the same compartment in Figure 12(e) is more or less uniform at 5,500 psi.

The effect of differential pressures is separately graphed in Figures 13(a)–13(d) with pressure profiles along four lines AB, CD, EF, and GH indicated in Figures 12(e) and 12(f). The pressure profiles along lines AB and CD (Figures 13(a) and 13(b)) are identical. These uniform pressure profiles are due to the symmetry of the initial conditions along the boundaries of the compartments in Figure 12(e). In contrast, for the case of Figure 12(f), pressure profiles across Compartments 2 and 3 vary laterally as follows from the two pressure profiles of Figures 13(c) and 13(d).

(a)

(b)

(c)

(d)

Figure 13: Pressure profiles along lines AB, CD, EF, and GH, where labels 1, 2, and 3 correspond to the pressure of the three compartments mentioned in Figure 12. The blue dots show the pressure value for each point along a specific line in the direction, as flow. (a) Pressure profile along line AB. (b) Pressure profile along line CD. (c) Pressure profile along line EF. (d) Pressure profile along line GH.

2.8.4. Case H: Multilayered Heterogeneous Porosity Zones with High Conductivity Gaps

Case H uses the same layered reservoir as for Case G but highlights the effect of heterogeneous porosity distribution on FP and TOF. The reservoir is still divided into layered zones, gray areas (Figure 12(a)) have higher porosity, 20%, and white areas have lower porosity, 5%. However, the permeability distribution across the layers is uniform with 100 mD.

The numerical streamline workflow is applied for both reservoir settings illustrated in Figures 12(a) and 12(b). The corresponding results of FP and TOF are plotted in Figures 14(a) and 14(b). Since the permeability distribution is uniform for both cases, the streamline trajectories are identical. Figure 14(a) shows that the flow runs faster through the two gaps with the decreased porosity and quickly arrives at the producer wells. In Figure 14(b), the second case with just one gap shows similar result for the two lower gaps, but the flow runs slower in the absence of the 2nd gap where no porosity change occurs.

(a)

(b)

Figure 14: Two central clay or shale walls obstruct flow, except for sand filled fractures that leave a high permeability conduit for the flood advance across the two central low permeability walls. Low permeability walls also occur above and below the well arrays. (a) The streamlines (red curves) and flood flight times (rainbow colors for TOF scaled in unit of days) for case of Figure 12(a), where sand porosity is 5%; clay is 20%; permeability everywhere equals 100 mD. (b) The streamlines (red curves) and flood flight times (rainbow colors for TOF scaled in unit of days) for case of Figure 12(b), where sand porosity is 5%; clay is 20%; permeability everywhere equals 100 mD.

3. Discussion

In this study, the relationships between the time of flight and flow path in porous media with spatial changes in permeability and porosity are derived and illustrated, both theoretically and numerically. Two fundamental rules are formulated to capture the essence of the relationships.

A principal conclusion of our study, based on careful derivations and systematically designed streamline simulations, is that permeability distribution alone governs the spatial distribution of streamlines. Porosity does have an influence on the time of flight but does not affect the streamline path. Those are not trivia conclusions, and although many practitioners may have empirical understanding of these relationships, we know of no systematic treatise with an explicit theoretical proof and substantiated numerical experimental validation for the following rules. Rule  1 states that an increase in permeability decreases the time of flight, whereas an increase in porosity increases the time of flight. Rule  2 in our present study states that the permeability uniquely controls the flight path of fluid flow in porous media; local porosity variations do not affect the streamline path.

While these two rules are generic and applicable to any porous medium, the motivation of our study mainly comes from the surging interest in shale oil and gas production, in particular when studying the effect of anisotropy fabrics on reservoir flow, to which our two rules apply. For example, the emergence of oil and gas recovery from shale formations has led to renewed interest in laboratory measurements of anisotropic permeability properties of shale samples from major shale plays (e.g., Eagle Ford, Barnett in the US, and emerging plays elsewhere, e.g., Longmaxi formation, Sichuan Basin, China, and Vaca Muerta, Argentina). Numerous studies have detailed the permeability properties (e.g., [37, 39]). For example, in Pan et al. [39], the authors measured that the shale reservoir could have a horizontal permeability (about 200 nD) that is 25 times the vertical permeability (about 8 nD), and Mokhtari and Tutuncu [37] showed that the horizontal permeability is almost 4 times the vertical permeability. More permeability measurement results are given in Table 2.

Our simulations of Cases E and F used permeability anisotropy ratios of 5 and 10, respectively, which are conservative anisotropy factors compared to the range of values given in Table 2. Our study highlights that permeability is indeed important for modeling of the fluid migration paths in any reservoir due to its impact on the spatial distribution of fluid velocities and henceforth control on the appearance of flow paths. Given the strong permeability anisotropy in shale layers all migration of hydrocarbons and other reservoir fluids through the rock matrix space will likely occur along flow paths aligned with the anisotropy plane as our experiments suggest that practically no flow can occur normal to such surfaces (Figure 11(b)). Our flow visualizations provide compelling evidence for our theoretical assertion that, in addition to the permeability, the porosity needs to be known accurately too, in order to be able to reliably model the time of flight in such reservoirs.

Our simulation of flow across sand compartments separated by shale walls (Figures 12 and 13) revealed that the flow conduits between such compartments act as pressure valves with steep pressure gradients and fast flow rate in conduits, whereas pressure in the sand compartments remain relatively uniform (Figure 12(e)). When more conduits connect to an adjacent high pressure zone upstream than to the lower pressure zone downstream, the central reservoir compartment will be partially overpressured (Figure 12(f)). This finding provides scope for future studies to model the occurrence of overpressured zones commonly encountered in shale-rich formations based on our streamline and pressure models of fluid migration in such zones. Such overpressured shales provide drilling hazards, both in onshore unconventional oil and gas plays [50–53] and in one of the principal reservoirs offshore (Wilcox formation, Gulf of Mexico) where clay barriers and baffles restrict fluid flow between sand compartments [54].

A related finding from the Appendix is that, assuming constant injector rates in the reservoir, the pressure increases with viscosity of the produced fluids. Pressure will also increase due to any decrease(s) in local permeability and/or porosity.

Further insight obtained from our flow simulations is that larger pore space increases residence time of fluid passing through each pore (slowing TOF in such high porosity zones). Our results also shed light on residence time of acid injectates. For example, acid stimulation will more effectively react with high porosity zones. A common conjecture is that a higher effect of acidization is due to the larger volume of acid or larger contact surface area available for the acid simulation. But with the result of our study a higher effect of acid could also be a consequence of longer residence time in the high porosity spaces. Also, acid injection may reach producer wells faster via the low porosity zones and the injected acid reaches the high porosity space much slower (Figure 8(c)). Additionally, such models need to account for the chemical reaction which affects the porosity and permeability and therefore could completely change the flow distribution over time.

The rules highlighted in our study also apply to streamline studies for well productivity optimization and history matching. Several studies have been completed by us applying streamline solutions to a variety of geothermal and hydrocarbon production systems, including the optimization of well locations and well production rates in order to obtain maximum water sweep across oil drainage volumes [13, 42, 55, 56]. For example, Weijermars et al. [13] studied the sweep efficiency of a line drive scheme for five injectors and five producers under different permeability and porosity distributions. In accordance with our current study, the flow will converge and run faster through some higher permeability zones, so that an optimization plan for injection rate was designed to even the water front for five different injectors to achieve the maximum sweep efficiency. Weijermars and van Harmelen [56] studied the advancement of sweep zones in waterflooding for reservoirs with impermeable faults. The effect of such faults is equivalent to a zone with infinitely low permeability, which will cause the streamline trajectories to detour around the fault.

4. Conclusions

High fidelity reservoir simulations that include complex reservoir geometry may obscure some of the systematic effects highlighted in our systematic study. Our conclusion is that permeability and porosity both need to be accounted for in well productivity models and in history matching process for the reservoir property parameters, for the following reasons:

(1)Spatial permeability changes affect the streamline pattern (or particle flight paths if flow is transient).(2)Spatial porosity changes do not affect the streamline pattern.(3)Porosity changes for otherwise unchanged permeability distributions will affect the time of flight, but not the streamline pattern.(4)When the permeability is anisotropic and/or heterogeneous, the streamline patterns are affected. However, any spatial porosity changes will affect the time of flight without changing the streamline patterns.(5)Any temporal permeability changes due to reservoir compaction as production proceeds will effectuate streamline migration over the life cycle of a field.

In addition to permeability and porosity gradients of otherwise isotropic domains, spatial anisotropy of permeability can independently affect well productivity. Anisotropy and permeability both control the flow (or flight) path, leading to dispersion or channeling (Figure 10) and asymmetry (Figure 11) of streamlines. Any spatial porosity jumps or gradients will affect the time of flight (Figures 6–9). Our synthetic flow visualizations may provide useful theoretical guidance for more advanced reservoir simulation studies, in which complex geometrical details and spatial patterns of heterogeneity and anisotropy could obscure and prevent recognition of some of the systematic effects illustrated in the present study.

Appendix

Relation between Pressure and Permeability/Viscosity

Here we demonstrate that when the fluid flux of the injectors and producers should remain constant, the reservoir pressure changes inversely proportional to the changes of permeability and proportionally to the changes in viscosity. We performed two sets of experiments. In the first set of experiments, we varied the permeability value from 100 mD to 50 mD and 25 mD. In the second set, the viscosity contrast of water and oil was kept constant but absolute values were changed from 0.5 cp/2 cp to 1 cp/4 cp and 2 cp/8 cp. The resulting pressure changes are summarized in Table 3; the corresponding pressure distributions are mapped in Figure 15. From Figure 15, we conclude that the pressure will change commensurately with changes of permeability and/or viscosity. Similar experiments have been completed to confirm that the same conclusion applies to the relation between pressure and porosity.

Table 3: Pressure changes according to different permeability and viscosity.

Porosity (%)
Permeability (mD)
Viscosity (water/oil) (cp)
Pressure distribution (psi)

0.2
100
0.5/2
362–7361

0.2
50
0.5/2
705–14590

0.2
25
0.5/2
1140–27691

0.2
100
1/4
705–14590

0.2
100
2/8
1140–27691

(a) 100 md; 0.5 cp/2 cp

(b) 50 md; 0.5 cp/2 cp

(c) 100 md; 1 cp/4 cp

(d) 25 md; 0.5 cp/2 cp

(e) 100 md; 2 cp/8 cp

Figure 15: The pressure distributions for different permeability and viscosity. The corresponding pressure range is given in Table 3. (a) Pressure distribution of the base case model with uniform permeability of 100 mD and the viscosity of water and oil is 0.5 cp and 2 cp, respectively. (b) Permeability is 50 md; other parameters are the same as in the original model. (c) Viscosity for water/oil is 1 cp/4 cp; other parameters are the same as in the original model. (d) Permeability is 25 md; other parameters are the same as in the original model. (e) Viscosity for water/oil is 2 cp/8 cp; other parameters are the same as in the original model.

Nomenclature

: | Porosity
: | Velocity vector
: | Velocity in the direction
: | Velocity in the direction
: | Velocity in the direction
: | Density
: | Pressure
: | Total compressibility
: | Dynamic viscosity
: | Permeability
: | Component of permeability tensor
: | Gravity in the direction
: | Volumetric flux in the direction
: | Initial gas flow rate
: | Stream function in 2D
: | Time of flight
: | Spatial distance along streamline.

SI Metric Conversion Factors

ft × 3.048e − 01 | =  m
  ×  2.832e− 02 | =  
cp × 1.0e − 03 | =  Pa·s
psi × 6.895e + 00 | =  kPa
mD × 9.869e − 16 | =  
nD × 1e − 6 | =  mD.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

Lihua Zuo gratefully acknowledges financial support from members of the Texas A&M University Joint Industry Project, MCERI (Model Calibration and Efficient Reservoir Imaging). This study was sponsored by start-up funds from the Texas A&M Engineering Experiment Station (TEES) provided to the corresponding author.

References

E. D. Pittman, “Relationship of porosity and permeability to various parameters derived from mercury injection-capillary pressure curves for sandstone,” AAPG Bulletin, vol. 76, pp. 191–198, 1992.

P. H. Nelson, “Permeability-porosity relationships in sedimentary rocks,” Log Analyst, vol. 35, no. 3, pp. 38–62, 1994.

A. Koponen, M. Kataja, and J. Timonen, “Permeability and effective porosity of porous media,” Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 56, no. 3, pp. 3319–3325, 1997.

A. Costa, “Permeability-porosity relationship: a reexamination of the Kozeny-Carman equation based on a fractal pore-space geometry assumption,” Geophysical Research Letters, vol. 33, no. 2, Article ID L02318, 2006.

P. C. Carman, “Fluid flow through a granular bed,” Transactions of the Institution of Chemical Engineers, vol. 15, pp. 150–156, 1938.

P. C. Carman, Flow of Gases through Porous Media, Butterworths Scientific Publications, London, UK, 1956.

J. Kozeny, “Uber kapillare leitung des wassers im boden,” Sitzungsberichte der Kaiserlichen Akademie der Wissenschaften, vol. 136, pp. 271–306, 1927.

G. I. Barenblatt, T. W. Patzek, V. M. Prostokishin, and D. B. Silin, “Oil deposits in diatomites: a new challenge for subterranean mechanics,” in Proceedings of the SPE/DOE Thirteenth Symposium on Improved Oil Recovery, pp. 899–907, Tulsa, Okla, USA, April 2002.

G. I. Barenblatt, M. Bertsch, and C. Nitsch, “Nonlocal damage accumulation and fluid flow in diatomites,” Communications in Applied Mathematics and Computational Science, vol. 1, pp. 143–168, 2006.

R. P. Batycky, A. C. Seto, and D. H. Fenwick, “Assisted History matching of a 1.4-million-cell simulation model for judy creek—a pool waterflood/HCMF Using a streamline-based workflow,” in Proceedings of the SPE Annual Technical Conference and Exhibition, Anaheim, Calif, USA, November 2007.

A. Datta-Gupta and M. K. King, Streamline simulation: Theory and Practice, SPE, Richardson, Texas, USA, 2007.

EDIAFILT, “Industrial utilization,” 2015, http://www.ediafilt.hu/en/teruletek.html.

R. Weijermars, A. van Harmelen, and L. Zuo, “Controlling flood displacement fronts using a parallel analytical streamline simulator,” Journal of Petroleum Science and Engineering, vol. 139, pp. 23–42, 2016.

J. M. Heidt and G. J. Follensbee, “Application of the Higgins-Leighton waterflood prediction technique to L-lalfway reservoirs in Northeastern British Columbia,” Journal of Canadian Petroleum Technology, vol. 10, no. 01, pp. 23–28, 1971.

J. C. Martin, P. T. Woo, and R. E. Wagner, “Failure of stream tube methods to predict waterflood performance of an isolated inverted five-spot at favorable mobility ratios,” Journal of Petroleum Technology, vol. 25, no. 02, pp. 151–153, 1973.

W. Hauber, “Prediction of waterflood performance for arbitrary well patterns and mobility ratios,” Journal of Petroleum Technology, vol. 16, no. 01, pp. 95–103, 1964.

R. V. Higgins, D. W. Boley, and A. J. Leighton, “Aids to forecasting the performance of water floods,” Journal of Petroleum Technology, vol. 16, no. 09, pp. 1076–1082, 1964.

J. L. LeBlanc and B. H. Caudle, “A streamline model for secondary recovery,” SPE Journal, vol. 11, no. 01, pp. 7–12, 1971.

J. Martin and R. Wegner, “Numerical solution of multiphase, two-dimensional incompressible flow using stream-tube relationships,” SPE Journal, vol. 19, no. 05, pp. 313–323, 1979.

J. H. Abou-Kassem and K. Aziz, “Analytical well models for reservoir simulation,” SPE Journal, vol. 25, no. 04, pp. 573–579, 1985.

D. O. Cox, “Waterflood performance estimation with a layered streamtube model,” in Proceedings of the Petroleum Industry Application of Microcomputers, Lake Conroe, Tex, USA, June 1987.

M. R. Thiele, Modeling multiphase flow in heterogeneous media using streamtubes [Ph.D. thesis], Stanford University, 1994.

M. R. Thiele, R. P. Batycky, M. J. Blunt, and F. M. Orr Jr., “Simulating flow in heterogeneous systems using streamtubes and streamlines,” SPE Reservoir Engineering, vol. 11, no. 01, pp. 5–12, 1996.

R. W. Parsons, “Directional permeability effects in developed and unconfined five-spots,” Journal of Petroleum Technology, vol. 24, no. 04, pp. 487–494, 1972.

RP. Batycky, MR. Thiele, RO. Baker, and SH. Chugh, “Revisiting reservoir flood-surveillance methods using streamlines,” SPE Reservoir Evaluation & Engineering, vol. 11, no. 2, pp. 387–394, 2008.

P. Samier, L. Quettier, and M. Thiele, “Applications of streamline simulations to reservoir studies,” in Proceedings of the SPE Reservoir Simulation Symposium, Houston, Tex, USA, February, 2001.

R. Raghavan and C. Chen, “Fractured-well performance under anomalous diffusion,” SPE Reservoir Evaluation & Engineering, vol. 16, no. 3, pp. 237–245, 2013.

R. Raghavan and C. Chen, “Rate decline, power laws, and subdiffusion in fractured rocks,” in Proceedings of the SPE Low Perm Symposium, Denver, Colo, USA, May 2016.

W. Rundell, X. Xu, and L. Zuo, “The determination of an unknown boundary condition in a fractional diffusion equation,” Applicable Analysis: An International Journal, vol. 92, no. 7, pp. 1511–1526, 2013.

Y. Luchko, W. Rundell, M. Yamamoto, and L. Zuo, “Uniqueness and reconstruction of an unknown semilinear term in a time-fractional reaction-diffusion equation,” Inverse Problems, vol. 29, no. 6, Article ID 065019, 2013.

Y. Luchko and L. Zuo, “-function method for a time-fractional reaction-diffusion equation,” Journal of Algerian Mathematical Society, vol. 1, pp. 1–15, 2014.

R. W. Holy and E. Ozkan, “A practical and rigorous approach for production data analysis in unconventional wells,” in Proceedings of the SPE Low Perm Symposium, Denver, Colo, USA, May 2016.

L. Zuo, J. Lim, R. Chen, and M. J. King, “Efficient Calculation of Flux Conservative Streamline Trajectories on Complex and Unstructured Grids,” in Proceedings of the 78th EAGE Conference & Exhibition, Vienna, Austria, 2016.

L. Zuo, W. Yu, and K. Wu, “A fractional decline curve analysis model for shale gas reservoirs,” International Journal of Coal Geology, vol. 163, pp. 140–148, 2016.

J. Bear, Dynamics of Fluids in Porous Media, Dover Publications, Mineola, NY, USA, 1988.

A. R. Bhandari, P. B. Flemings, P. J. Polito, M. B. Cronin, and S. L. Bryant, “Anisotropy and stress dependence of permeability in the Barnett Shale,” Transport in Porous Media, vol. 108, no. 2, pp. 393–411, 2015.

M. Mokhtari and A. N. Tutuncu, “Characterization of anisotropy in the permeability of organic-rich shales,” Journal of Petroleum Science and Engineering, vol. 133, pp. 496–506, 2015.

O. Kwon, A. K. Kronenberg, A. F. Gangi, B. Johnson, and B. E. Herbert, “Permeability of illite-bearing shale: 1. Anisotropy and effects of clay content and loading,” Journal of Geophysical Research: Solid Earth, vol. 109, no. B10, Article ID B10205, 2004.

Z. Pan, Y. Ma, L. D. Connell, D. I. Down, and M. Camilleri, “Measuring anisotropic permeability using a cubic shale sample in a triaxial cell,” Journal of Natural Gas Science and Engineering, vol. 26, pp. 336–344, 2015.

P. J. Armitage, D. R. Faulkner, R. H. Worden, A. C. Aplin, A. R. Butcher, and J. Iliffe, “Experimental measurement of, and controls on, permeability and permeability anisotropy of caprocks from the CO₂ storage project at the Krechba Field, Algeria,” Journal of Geophysical Research: Solid Earth, vol. 116, no. b12, Article ID B12208, 2011.

D. W. Pollock, “Semianalytical computation of path lines for finite‐difference models,” Groundwater, vol. 26, no. 6, pp. 743–750, 1988.

R. Nelson, L. H. Zuo, R. Weijermars, and D. Crowdy, “Outer boundary effects in a petroleum reservoir (Quitman field, east Texas): applying improved analytical methods for modelling and visualization of flood displacement fronts,” Journal of Petroleum Science and Engineering, moderate revision, 2017.

R. Weijermars, A. van Harmelen, and L. Zuo, “Flow interference between frac clusters (Part 2): field example from the Midland basin (Wolfcamp formation, Spraberry trend field) with implications for hydraulic fracture design,” in Proceedings of the Model Verification and Basic Cases, Presented at the Unconventional Resources Technology Conference, SPE-2670073B, pp. 989–1004, Austin, Tex, USA, 2017.

R. Weijermars, A. van Harmelen, L. Zuo, I. N. Alves, and W. Yu, “High-resolution visualization of flow interference between frac clusters (part 1): model verification and basic cases,” in Proceedings of the Model Verification and Basic Cases, Presented at the Unconventional Resources Technology Conference, SPE-2670073A, pp. 963–988, Austin, Tex, USA, July 2017.

R. Weijermars, L. Zuo, and I. Warren, “Modeling Reservoir circulation and economic performance of the Neal Hot Springs geothermal power plant (Oregon, U.S.): an integrated case study,” Geothermics, vol. 70, pp. 155–172, 2017.

M. Muskat, Physical Principles of Oil Production, McGraw-Hill, New York, NY, USA, 1949.

M. Muskat, “The theory of potentiometric models,” Transactions of the AIME, vol. 179, no. 01, pp. 216–221, 1949.

F. Trippetta, B. M. Carpenter, S. Mollo, M. M. Scuderi, P. Scarlato, and C. Collettini, “Physical and transport property variations within carbonate-bearing fault zones: insights from the monte maggio fault (Central Italy),” Geochemistry, Geophysics, Geosystems, pp. 1–16, 2017.

J. Bear, Dynamics of Fluids in Porous Media, Elsevier, New York, NY, USA, 1972.

R. A. Lane and L. A. Macpherson, “A review of geopressure evaluation from well logs—Louisiana Gulf Coast,” Journal of Petroleum Technology, vol. 28, no. 09, pp. 963–971, 1976.

M.-K. Lee and D. D. Williams, “Paleohydrology of the Delaware Basin, Western Texas: overpressure development, hydrocarbon migration, and ore genesis,” AAPG Bulletin, vol. 84, no. 7, pp. 961–974, 2000.

M. Luo and M. Baker, “Distribution and generation of overpressure system, eastern delaware basin, Western Texas, and Southern New Mexico: reply,” AAPG Bulletin, vol. 79, no. 12, pp. 1823-1824, 1995.

M. Luo, M. Baker, and D. LeMone, “Distribution and generation of overpressure system, eastern delaware basin, Western Texas, and Southern New Mexico,” Journal of Petroleum Technology, vol. 78, no. 9, pp. 1386–1405, 1994.

H. M. Li, C. Genty, T. Sun, and J. K. Miller, “Modeling flow barriers and baffles in distributary systems using a numerical analog from process-based model,” AAPG: Search and Discovery, Article ID 50164, 2009.

R. Weijermars and A. van Harmelen, “Breakdown of doublet recirculation and direct line drives by far-field flow in reservoirs: implications for geothermal and hydrocarbon well placement,” Geophysical Journal International, vol. 206, no. 1, pp. 19–47, 2016.

R. Weijermars and A. van Harmelen, “Advancement of sweep zones in waterflooding: conceptual insight based on flow visualizations of oil-withdrawal contours and waterflood time-of-flight contours using complex potentials,” Journal of Petroleum Exploration and Production Technology, vol. 7, no. 3, pp. 785–812, 2017.