This blog post describes a new capability in ResFrac to capture the effect of ‘fracture swarms’ on production decline trends. Based on work from Acuna (2020), the idea is that variable spacing between fractures causes a gradual onset of production interference. Fractures in a swarm may be numerous and tightly spaced,  so rather than representing each individual crack in the model, we treat each swarm as a single crack and use a numerical technique to capture their effects. In ResFrac, this capability is useful because it provides another mechanism for explaining (and matching) production drawdown trends. For further details, refer to Section 19.10 from McClure et al. (2022).
 

Background – Deviation From Linear Flow

‘Linear flow’ is a flow regime in which pressure drawdown is primarily driven by fluid flowing directly inwards towards the propped hydraulic fracture(s), such that the flow geometry is approximately perpendicular to the fracture face, or ‘linear’ (shown in the image below). During production in a multifractured shale well, the fractures are long and the permeability is low, so ‘linear flow’ is a dominant flow geometry.

 
With linear flow, if we produce at constant pressure, then rate should decay with the square root of time. Under these conditions, if we plot the reciprocal of production rate against the square root of time, we will observe a straight line.

In real production data, the bottom-hole pressure changes over time, so we cannot simply plot the reciprocal of rate versus time. Instead, a specialized plot of ‘reciprocal productivity index’ (RPI) versus ‘square root of material balance time’ can be constructed (Liang et al., 2011). This plot shows approximately what the data ‘would have’ looked like if the production had occurred at constant BHP.

Very often, in real data, we observe a loss of production that scales more rapidly than would be expected from pure linear flow.

The figure below (source) shows examples with (a) pure linear flow, and (b) deviation from linear flow. Linear flow corresponds to a straight line. When productivity falls off more rapidly than expected from pure linear flow, the plot bends upwards (deviation from linear flow).

 
What causes deviation from linear flow? The ‘classical’ cause is fracture-to-fracture interference, which happens when the drainage volumes from adjacent fractures begin to overlap. If we attempt production history matching and assume a single dominant fracture per cluster, then deviation from linear flow often occurs far too early to be caused solely by fracture-to-fracture interference. A variety of processes have been proposed as potential causes: multiphase flow, pressure-dependent system permeability, pressure/stress-dependent fracture conductivity, time-dependent fracture conductivity loss, and fracture swarming (Clarkson et al., 2013; Acuna, 2020; Carlsen et al., 2021; Carlsen and Whitson, 2022).

It is often observed that the productivity decline scales according to an equation:

Where RPI is reciprocal productivity index (reservoir pressure minus BHP divided by production rate), t is time, and α is a coefficient. Linear flow corresponds to α of 0.5. Often it is observed that productivity decays with an α value greater than 0.5, corresponding to a more-rapid loss of productivity (ie, deviation from linear flow).

In ResFrac (or any numerical simulator), it is relatively straightforward to include multiphase flow effects, time-dependent conductivity loss, and pressure/stress-dependent permeability and conductivity. However, the potential effect of ‘fracture swarming’ involves sub-model-scale processes, and there is no standard way of incorporating it into a simulator.
 
 

Incorporating ‘fractal fracture swarms’ into ResFrac

Acuna (2020) proposes that fracture swarming is a leading cause of deviation from linear flow (α greater than 0.5). Observations from core through studies have shown that fractures often propagate in swarms, rather than as single individual cracks (Warpinski et al., 1993; Raterman et al., 2019; Gale et al., 2018; Gale et al., 2022). The practical implications of these ‘fracture swarm’ observations are still being debated.

Raterman et al. (2019) suggests that only a fraction of the total hydraulic fracture strands observed in the core contain proppant and contribute to the long-term drainage of fluid towards the well. Also, the Raterman et al. (2019) observations suggest that these major draining fractures are spaced erratically. Sometimes, there may be only a few feet between a draining fracture, and in other cases, there may be tens of feet of gap between flowing fractures.

When the drainage regions from adjacent fractures begin to overlap, this causes a deviation from linear flow (an upward bend in the RTA plot). If there was a single fracture per perforation cluster and they were uniformly spaced, then the drainage regions from the adjacent fractures would all interfere at roughly the same time, because they would all be the same distance apart. This would result in a relatively long ‘linear flow’ period followed by a relatively abrupt upward bend on the RTA plot. However, if the fractures are spaced erratically, they will interfere at different times, causing an earlier deviation from linear flow, and a more gradual onset. Also, if the fractures form in swarms, then there may be an early time ‘flush production’ as closely spaced fractures rapidly drain the rock between them.

Acuna (2020) points out that the effect of variable fracture spacing can be modeled as a ‘fractal fracture swarm.’ The ‘fractal dimension,’ D, quantifies the degree of clustering and/or uniformity in the fracture spacing. A particular value of D implies a particular value of α in the rate-transient (2016; 2018; 2020). Acuna (2020) validates his derivations with direct numerical solutions of arrays of cracks with spacing described by different ‘fractal dimensions.’

If a ‘single crack’ is actually representing a ‘swarm of fractures,’ this will increase the effective fracture surface area. If there are ten cracks in the swarm, the initial effective fracture area will be 10x higher. Hypothetically, let’s assume that some of those cracks are within only a few inches of each other. Once the drainage volumes of these fractures have overlapped (which will happen relatively soon), they will subsequently drain from outside their combined area.

Further, let’s assume that all of these fractures are within a region 10 ft wide. Once the drainage volumes of all the fractures within the 10 ft band have merged, they will now be draining only from the rock outside the 10 ft band, and their effective drainage volume at distances greater than 10 ft will now be the same as if they were only a single crack.

As described by the ‘butterfly plot’ from Acuna (2018), the net effect of fracture swarming is that the effective drainage volume decreases with distance from the fractures. Because the drained volume spreads out from the fracture(s) over time, this means that the effective draining area decreases with both distance and time, which causes a reduction in productivity and an upward deviation from linear flow on the RTA plot.

How should this process be modeled in a numerical simulator? One option would be to individually represent and mesh every individual crack. However, if the fracture swarms involve a large number of closely spaced fractures, this would require a very detailed discretization, at high computational cost.

Instead, we can treat each ‘swarm’ as a single crack, but capture the effect of fracture swarming by incorporating calculations at the submodel scale. During fracture propagation, this can be accomplished by scaling toughness, viscous pressure drop, and leakoff (Fu et al., 2020; Section 17 from McClure et al., 2022). For production, the ‘fractal fracture swarm’ concept from Acuna (2020) is suitable for capturing these effects. However, Acuna (2020) uses analytical solutions; we need a way of implementing this concept into ResFrac, which is a general numerical simulator. With testing, I’ve found that this can be accomplished by integrating with ResFrac’s ‘1D submesh method’ (McClure, 2017).

The ‘1D submesh’ approach is used in ResFrac to accurately capture fluid exchange between matrix and fracture elements, while maintaining computational efficiency. For every ‘connection’ between a fracture element and a matrix element, a function of ‘rock volume versus distance’ is constructed. To construct these functions, each matrix element (containing fracture elements) is subdivided into smaller volumes, and each is assigned uniquely to one of the submeshes corresponding to a fracture-matrix connection. The volumes assigned to each submesh are then used to generate the ‘volume versus distance’ curves used by the submesh.

The approach transforms arbitrary flow geometries into a generalized one-dimensional problem. For example, with ‘radial flow,’ the volume of rock within a certain distance r of the center scales with the perimeter (i.e., scales linearly with r), and the total volume within a certain distance scales with the integral, or r². With ‘linear flow,’ the volume of rock at a certain distance, x, from the fracture is constant, and the cumulative volume within a certain distance scales linearly.

Next, let’s consider the case of production into a vertical, hydraulically fractured well. In this case, as is known from classical PTA/RTA, flow is initially linear, then elliptical, and eventually radial. So, in this case, a function of ‘volume of rock’ versus ‘flow distance’ will initially be constant (corresponding to linear flow), but then transition into scaling with distance (corresponding to radial flow).

The ‘1D submesh’ approach is related to methods such as streamline simulation, the fast-marching method (Datta-Gupta et al., 2011), and the butterfly model described by Acuna (2018).

Once the submesh function has been constructed, it can be meshed with very fine mesh refinement near the fracture, which enables solution accuracy even at short timescales with very low matrix permeability. 1D problems lead to banded diagonal systems, which can be solved very efficiently.

The 1D submesh approach is advantageous because: (a) it avoids the need for a conforming mesh between the fracture and matrix elements, which would require complex and frequent unstructured remeshing during propagation, (b) it maintains accuracy at all timescales, even with very low matrix permeability, and (c) it adds only modestly to the overall computational cost.

In contrast, the ‘embedded discrete fracture method’ (EDFM) uses a nonconforming mesh to handle fracture-matrix fluid flow (Li and Lee, 2008), but does not use a ‘submesh’ to capture small-scale effects. As a result, it can suffer large numerical errors at relatively short timescales and/or with very low matrix permeability.

The 1D submesh very naturally incorporates the effect of fractal fracture dimension. Once the 1D submesh has been constructed (a function describing drainage volume versus distance for each fracture-to-matrix connection), we then adjust the volumes as a function of distance to account for fracture swarming.

ResFrac permits the user to specify a ‘fractal dimension’ (referred to as D; following Acuna 2020). If you set to 0, it has no effect. If you set to 0.3, there is a mild effect. If set to 0.5-0.7, the effect is strong.

The setting is implemented by increasing the flow ‘volume’ and ‘area’ in the vicinity of the fracture, with the effect diminishing with distance, until 20 ft from the fracture face, at which point there is no longer any adjustment.

Let’s say, for example, that the baseline 1D submesh is defined by a function of distance from the fracture, A(x), that defines the cross-sectional area for flow as a function of distance (or equivalently, the incremental volume per distance). If you specify ‘submesh fractal D,’ then the function is modified to a new function, Ã(x), according to the equations:

The length of the submesh may be truncated to ensure that it does not ‘overcount’ the total volume of the element. For example, let’s say that the total volume of the original, unmodified submesh is:

where L is the length of the submesh from the fracture face. The modified submesh must have the same volume, and so the submesh is truncated at a distance $\tilde{L}$, such that:

The figures below show the effect of using ‘submesh fractal D,’ using values of 0, 0.3, and 0.6. The ‘submesh fractal D’ setting increases early-time production. It is equivalent to assuming that each individual ‘crack’ represented in the simulator is actually a swarm of fractures. With higher ‘fractal D,’ each swarm has a larger number of fractures. As a result, simulations with higher value have more production. However, with more fractures, interference between fractures occurs more rapidly, which is why there is more curvature on the RTA plot.

 

 

 

The fractal swarm calculation increases early-time production, and so when increasing D, it is often useful to simultaneously decrease system permeability.

When matching the RTA curvature in real data, you have a variety of tools – relative permeability loss with drawdown, pressure dependent permeability loss, pressure/stress dependent conductivity loss, and time-dependent conductivity loss. The ‘fractal D’ setting gives one more tool that can be used to match RTA curvature. It may often be useful (and realistic) to include multiple simultaneous processes. For example, the figure below shows the same simulation as shown above with ‘fractal D’ set to 0.6, but also with some time-dependent conductivity loss. The simulation exhibits even stronger RTA curvature.

 
To validate numerical accuracy, I ran linear-flow simulations with constant production rate and different values of D. Then, I made log-log plots of t*dP/dt versus time. The slope of the line on the log-log plot is equal to the value of α. According to the derivations from Acuna (2020), for linear flow geometry, the slope α observed on the plot should be equal to 1-δ, where δ=(1-D)/2.

With D equal to zero, the log-log plot should have a slope of 0.5, corresponding to standard linear flow:

 

With D of 0.2, the slope should be 0.6:

 

With D of 0.4, the slope should be 0.7:

 
With D of 0.6, the slope should be 0.8:

 

To estimate the slopes of the t*dP/dt curves (shown in purple in the upper right panels), pick two points that are separated by one log-10 cycle. The log slope is equal to the log-10 of their ratio. The plots above show the expected slopes, given the specified values of D.

The lower right plots show that the matrix region in this simulation meshed with just one element in the direction perpendicular to the fractures. The visualization only shows the ‘average’ pressure in the matrix elements. There isn’t a capability to visualize the details of the pressure drawdown ‘within’ the 1D submesh in those elements. Thus, in the actual calculation used by the simulator, the 1D submesh grid incorporates the small-scale details of the pressure distribution, even though this isn’t apparent in the visualization.
 
 

Recap

A variety of processes can cause deviation from linear flow during shale production. Fracture swarms have been proposed as a physically plausible mechanism that could cause deviation from linear flow behavior seen in real data (Acuna 2020).

We’ve implemented the ‘fracture swarm’ effect in ResFrac with a ‘fractal dimension’ adjustment to our ‘1D submesh method.’ The result is a useful tool in the toolbox for history matching and, specifically, for matching long-term RTA trends.
 
 

References

Acuna, Jorge A. 2016. Analytical Pressure and Rate Transient Models for Analysis of Complex Fracture Networks in Tight Reservoirs. Paper URTeC 2429710 presented at the Unconventional Resources Technology Conference, San Antonio, TX.

Acuna, Jorge A. 2018. Straightforward Representative Fluid Flow Models for Complex Fracture Networks in Unconventional Reservoirs. Paper URTec-2876208 presented at the Unconventional Resources Technology Conference, Houston, TX.

Acuna, Jorge. 2020. Rate transient analysis of fracture swarm fractal networks. Paper URTeC-2020-2118 presented at the Unconventional Resources Technology Conference, Austin, TX.

Carlsen, Mathias Lia, Braden Bowie, Mohamad Majzoub Dahouk, Stian Mydland, Curtis Hays Whitson, and Ilina Yusra. 2021. Numerical RTA in tight unconventionals. Paper SPE-205884-MS presented at the Annual Technical Conference and Exhibition, Dubai, UAE.

Carlsen, Mathias Lia, and Curtis Hays Whitson. 2022. Numerical RTA extended to complex fracture systems: Part 2. Paper SPE-210420-MS presented at the Annual Technical Conference and Exhibition, Houston, TX.

Clarkson, C. R., F. Qanbari, M. Nobakht, and L. Heffner. 2013. Incorporating geomechanical and dynamic hydraulic-fracture-property changes into rate-transient analysis: Example from the Haynesville shale. SPE Reservoir Evaluation & Engineering 16(3), 303-316.

Datta-Gupta, Akhil, Jiang Xie, Neha Gupta, Michael J. King, and W. John Lee. 2011. Radius of investigation and its generalization to unconventional reservoirs. Journal of Petroleum Technology.

Fu, Wei, Joseph Morris, Pengcheng Fu, Jixiang Huang, Chris Sherman, Randolph Settgast, Hui Wu, Frederick Ryerson. 2020. Developing upscaling approach for swarming hydraulic fractures observed at Hydraulic Fracturing Test Site through Multiscale Simulations. Paper SPE-199689-MS presented at the SPE Hydraulic Fracturing Technology Conference and Exhibition, The Woodlands, TX.

Gale, Julia F. W., Sara J. Elliott, Stephen E. Laubach. 2018. Hydraulic fractures in core from stimulated reservoirs: Core fracture description of the HFTS slant core, Midland Basin, West Texas. Paper URTeC: 2902624.

Gale, Julia Fiona Wells, Sara J. Elliott, Bethany G. Rysak, and Stephen E. Laubach. 2022. The critical role of core in understanding hydraulic fracturing. In: Neal, A., Ashton, M., Williams, L.S., Dee, S.J., Dodd, T.J.H. and Marshall, J.D. (eds) Core Values: the Role of Core in Twenty-first Century Reservoir Characterization. Geological Society, London, Special Publications, 527.

Liang, P., L. Mattar, and S. Moghadam. 2011. Analyzing variable rate/pressure data in transient linear flow in unconventional gas reservoirs. Paper CSUG/SPE 149472 presented at the Canadian Unconventional Resources Conference, Calgary, Alberta, Canada.

McClure, Mark. 2017. An accurate and efficient method for calculating fluid exchange between fractures and matrix with a non-confirming mesh. arXiv:1709.02493.

McClure, M., C. Kang, S. Medam, C. Hewson, and E. Dontsov. 2022 ResFrac Technical Writeup. [Online]. Available: https://arxiv.org/abs/1804.02092.

Li, L., and S. H. Lee. 2008. Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media. SPE Reservoir Evaluation and Engineering 11(4): 750-759, doi:10.2118/103901-PA.

Raterman, Kevin T., Yongshe Liu, and Logan Warren. 2019. Analysis of a drained rock volume: An Eagle Ford example. Paper URTeC-2019-263 presented at the Unconventional Resources Technology Conference, Denver, CO.

Warpinski, N. R., J. C. Lorenz, P. T. Branagan, F. R. Myal, and B. L. Gall. 1993. Examination of a cored hydraulic fracture in a deep gas well. SPE Production & Facilities 8(3), 150-158.