Mark McClure and Ankush Singh

Introduction

In this blog post, we review a refracturing and economic optimization case study. The model and history match are loosely based on “SPE Data Repository Well #1,” a publicly available refracturing and production dataset from the Eagle Ford shale. This post extends the analysis that we presented at a recent SPE workshop, “What New for PTA and RTA”.

The initial completion used 300 ft stage length with 100 ft cluster spacing. The refrac was performed by running a new liner and reperforating at tighter spacing. The refrac completion used 210 ft stage length with roughly 21 ft cluster spacing. The initial job was pumped with about 800 lbs/ft of proppant and 16 bbl/ft of fluid. The refrac was performed with roughly 2900 lbs/ft of proppant and 30 bbl/ft of fluid.

When history matching post-refrac production, it was critical to calibrate parameters related to the amount of crossflow ‘outside casing’ from the new perforations back into the original fractures. This was necessary because if the model assumed that the new fractures were entirely isolated from the original fractures, it overpredicted the actual production post-refrac.

In the economic optimization, the optimal design was strongly dependent on the choice of ‘objective function’. Maximization of NPV/section resulted in a much more aggressive well spacing than maximization of discounted return on investment (DROI). We will expand on this topic in an upcoming blog post, discussing how frac design and well spacing can be connected to economic optimization of overall field development.

The analysis in this blog post was performed with more limited data than would typically be available in a modeling project. For many key inputs, values were selected to be ‘reasonable,’ and were not taken directly from the actual dataset. In some cases, the data was incompletely provided, and in a few cases, we used different values from the information specified in the original dataset. Furthermore, the economic deck used in the optimization was entirely hypothetical. Because of these simplifications, the specific results in this case study should be not considered ‘realistic.’ This post is only a demonstration of the workflow. It is not a concrete, real-life optimization of an actual frac design. A full analysis based on more detailed data and realistic economics numbers may have arrived at significantly different results.

 

Model setup

The simulation includes 630 ft along the lateral. The results are scaled proportionally for full-well production. The actual lateral length was 5778 ft. This ‘actual’ well-length was used in the history match; however, in the economic optimization calculations, we assumed 10,000 ft lateral length. In all cases, the simulations were run out to 30 years.

In the history match, it was assumed that there are not any nearby offset wells. When modeling the refrac in the economic optimization, it was assumed that infill wells were drilled at the same time as the refrac. Thus, in the economic optimization, ‘well spacing’ was a key optimization variable.

We used formation properties provided with the SPE Data Repository Well #1. For properties that were not specified, we assumed what we felt were ‘reasonable’ values for the Eagle Ford shale.

The model has a significant stress barrier above the Eagle Ford shale, confining fractures within the zone. The fractures extend to a maximum height of about 400 ft. The well is landed in an 80 ft thick ‘main pay’ zone with higher permeability than the surrounding formation. In the ‘main pay’ zone, the effective permeability to oil in the history matched model (including the effect of rel perm) is 611 nd. In the surrounding formation, the effective permeability to oil is 153 nd. The image below shows the stress and permeability layering.

The porosity is set to 6%, and the oil saturation is set to 0.74. Young’s modulus varies from 3.5e6 to 4.5e6 psi. The relative permeability curves are simple Brooks-Corey curves, with exponents of 2 for water and oil and 3 for gas. The residual gas saturation is set to 0.03. The residual water and oil saturations were varied as part of the history matching process, and ultimately set to 0.44 and 0.15. The pore pressure gradient is set to 0.44 psi/ft, yielding a formation pressure of 3400 psi (lower than the pressure that was actually specified in the SPE Data Repo Well #1 dataset).

The initial producing GOR is 336 scf/STB, and the bubble point pressure is 1211 psi. We used the ResFrac black oil wizard to generate a black oil table from these values, based on correlations.

The history match set the residual water saturation higher than the initial water saturation, meaning that there was no mobile water in the system at initial conditions. This is consistent with the production data – even after a few years of production, the actual produced water volume was only about 25% of the total injected volume. Flowback of injected fluid from the formation, after it had leaked off, was modeled using the ‘water banking’ capability within ResFrac.

The system permeability irreversibly decreases with depletion, following an exponential curve, with up to a 300x decrease with drawdown to the minimum BHP of 600 psi. In addition, the model assumes both effective stress and time-dependent conductivity loss in the proppant pack, following the relations described in the ResFrac technical writeup (McClure et al., 2022). After two months of depletion, the proppant pack conductivities are on the order of 25 md-ft. After 2 years of depletion, the fracture conductivity has decreased by roughly a factor of ten, to around 2.5 md-ft. The time-dependent conductivity causes conductivity to drop even further over the long-term. Ten years after the refrac, the propped conductivity is on the order of 0.025 md-ft. After 30 years, propped conductivity is on the order of 0.001 md-ft.

The actual production data is shown in the figure below. On the y-axes, the production volumes correspond to the 630 ft length of the sector model, not the full-well production. The right panel shows an RTA plot of the production results prior to the refrac.

In the actual data, the first two months of production were performed with BHP around 2600 psi. Subsequently, BHP dropped to 630 psi (presumably due to the installation of an ESP), and then remained roughly flat. The simulations were performed by specifying BHP and simulating rate.

In the RTA plot, the first two months of production (which had high BHP) are apparent in the early, low-steepness straight line. When BHP drops with the installation of the BHP, production rate does increase. However, the RTA curve shifts up to a much steeper trend, suggesting either lower effective fracture area or lower permeability. This phenomenon was matched in the model with the stress and pressure-dependent conductivity and system permeability. The match to the long-term upward curvature of the RTA trend was assisted by the time-dependent conductivity loss mechanism.

 

History matching

The automated history matching algorithm was set up to vary five parameters: residual water saturation, residual oil saturation, a global permeability multiplier, and two parameters related to crossflow between fractures during the refrac (discussed below). We would ideally have data constraining fracture length and height (such as from microseismic, sealed-wellbore pressure monitoring, and/or offset fiber). Also, we would ideally have interference tests to help us understand effective producing/propped length. In the absence of this data, for this ‘hypothetical’ case study, we chose parameters that would yield reasonable values for the Eagle Ford, based on published results, such as from Brinkley et al. (2021).

In the ResFrac automated history matching tool, you can set up ‘parameter groups’ that are varied to achieve the match. A parameter group may include just one parameter (such as the global permeability multiplier) or a group of parameters (such as the permeability of individual layers). For each parameter group, you specify min/max ranges for multipliers (linear or logarithmically spaced) or adders, which are applied uniformly to each parameter in the group.

The algorithm was set up to match three objective functions: the time-series of cumulative oil versus time, the time-series of cumulative water versus time, and the cumulative oil at one particular point in time – at the end of production prior to the refrac. Each of the three objectives was weighted equally. The third objective was partly redundant with the first. This was intentional; it used to emphasize matching cumulative oil production over the six months prior to the refrac. Often, we would include a time-smoothed GOR trend in the history match. However, in this particular dataset, GOR was quite low, and it did not have a very distinctive or interesting trend. Thus, for simplicity, it was left out of the calibration.

 

History matching the refrac production data

It was fairly straightforward to match the initial production data. However, matching the refrac production required some relatively ‘advanced’ features. Prior to the refrac, a new linear was run, cementing off the original perf clusters, and then new perforations were shot at a tighter spacing. This is handled in ResFrac by specifying ‘deactivation’ times and ‘activation’ times for different clusters. Here’s the challenge – if we assume that a new fracture forms at each new perf cluster, and that they are perfectly isolated from the fractures created during the original frac job, then the model will tend to overpredicts production post-refrac.

It is widely recognized that flow outside casing (either through the cement sheath or an axial fracture) can cause crossflow between fractures. Further, because of the passage of time, the cement outside the original casing is quite likely to be degraded. There is also a possibility that the newly formed fractures bend towards the original fracs due to stress rotation (Rezaei et al., 2017). In the past, I’ve written about how I am skeptical that there is large scale ‘turning’ of hydraulic fractures due to depletion. I am skeptical because: (a) there is a lack of observational data to support the hypothesis and (b) analyses predicting turning usually rely on greatly underestimated values for horizontal stress anisotropy. However, a refrac is a special case – the new fractures initiate very close to the original depleted fractures. In a refrac, it is plausible to suspect that some of these fractures may tilt over and collide with the original fractures.

These crossflow processes are modeled explicitly in ResFrac. The user specifies a ‘cased well and fracture collision distance.’ If a fracture physically intersects a well within that distance of a perforation, it forms a hydraulic connection directly into that perforation. These connections occur both for fractures propagating from other wells and for fractures initiating from the same well.

In the simulations for this demo, we set the ‘cased well and fracture collision distance’ to 80 ft. This is less than the original cluster spacing (100 ft), but sufficiently wide that there is significant crossflow from the ‘refrac’ perfs into the original fractures, and between each other (at 21 ft spacing). Limited-entry completion causes flow out of the well during the refrac to be fairly uniform. However, once outside the casing, the fluid is able to crossflow between the clusters. This defeats limited-entry and allows fluid to localize into a smaller number of fractures, as if the job had been pumped without  limited entry.

If we run the refrac simulation with completely unrestricted crossflow between the new perfs and the original fractures (and each other), we find that the simulation underpredicts  production because there is too much  localization of flow. Thus, we need to have some kind of intermediate case – where crossflow can occur, but it is not entirely unrestricted.

ResFrac controls the strength of crossflow between clusters with a handful of parameters. First, the parameter “Cased well and fracture connection additional NW deltaP” affects the ability of fluid to flow from the perforations into other fractures via the crossflow mechanism. For such connections, the code introduces a ‘near-wellbore tortuosity’ flow term that is proportional to the square root of flow rate. Optionally, this ‘additional NW deltaP’ coefficient can be made a function of distance. However, this was not done in the simulations for this case study.

This ‘additional NW deltaP’ parameter defines the coefficient for the near-wellbore pressure drop term. Second, the code permits crossflow directly between the fracture elements connected to the well. These crossflow terms are restricted by a ‘flow barrier’ treated like flow in a series (harmonic averaging) and with strength proportional to the distance between the fractures. The strength of the flow barrier is given by the parameter “Connect frac through ‘cased well fracture collision distance’ transmissibility barrier,” which is given in units of md-ft.

So, to conclude, the strength of crossflow between the fractures is greater when the ‘additional NW deltaP’ term is lower and the ‘transmissibility barrier’ term is higher. With zero crossflow, the model production is too high, and that with unlimited crossflow, the model production is too low. Thus, we need to find the values of the crossflow parameters that are intermediate between the extremes and that provide a match to the data. This is accomplished by the automated history matching algorithm.

History matching results

The history matching algorithm runs three generations of simulations. First, it runs a set of simulations that broadly span the range of possible parameters. Then, it fits a fast-running proxy function to the results. From the proxy model, it creates a second set of simulations – based both on places where the proxy model has uncertainty (to make sure it doesn’t miss anything) and where the results look most promising. After the second set of simulations, the process is repeated a third time, and finally, the results are ranked from best to worse. Using the final proxy model, we can derive marginal distributions for parameters, to assess nonuniqueness (Kang et al., 2002).

The figure below shows the highest-ranked result from the history matching algorithm. The dashed lines are the ‘actual’ data and the solid lines are the simulated results.

The figure below shows the fracture geometry, proppant distribution, drawdown, and production distribution along the fractures after six years of simulation time.

For comparison, here is the same figure right before the refrac:

We can see that the refrac created a significant number of new fractures. It also pumped a significant amount of fluid back into the original fractures. Much of the fluid pumped into the original fractures was ‘inefficient’ because it added propped to already-depleted rock. However, because the refrac had so much more proppant than the original frac job, it actually extended the propped length away from those fractures and accessed new rock.

The maximum propped half-length was 1200 ft. It should be noted – in many projects, we have interference tests that help us constrain propped half-length. Propped length is affected by the parameter “maximum proppant immobilization,” which captures the tendency of proppant to be held up and trapped as it flows through a fracture (Maity and Ciezobka, 2020). Larger amount of proppant lead to shorter propped fracture lengths. In the absence of data to calibrated propped length, we used 0.2 lbs/ft^2 for “maximum proppant immobilization,” which is an ‘intermediate’ value. Some formations exhibit lower trapping, and others exhibit significantly higher.

The table below summarizes the final results from the automated history match.

 

Economic optimization

With the history-matched model, we can now pursue economic optimization. Let’s assume that nearby, we have other, similar parent wells that we would like to refrac. Based on the data from the history matching dataset, what would be the best design? Further, let’s assume that at the same time as the refrac, infill wells are drilled on either side, such that the well experiences production interference from neighboring wells after the refrac.

We would like to know – what is the optimal job size and well spacing around the parent well?

The effect of ‘well spacing’ is mimicked with a ‘zero perm cube.’ When the well is put on production, no flow boundary conditions are placed on either side, at a distance equal to half the well spacing. The no flow boundary condition prevents ‘double counting’ of production that could have occurred from neighboring wells.

We could have chosen to include more than one well in the optimization, in order to account for offset wells. However, the ‘zero perm cube’ approach works reasonably well, and so for simplicity, this was the approach used in this case study.

The optimization algorithm works similarly to the history matching algorithm. It runs a series of simulations, fits a proxy model, and then repeats the process through two or three iterations.

The optimization assumed: (a) \$70/bbl oil; (b) initial drilling and ‘fixed’ completion cost of \$4M; (c) marginal fracturing cost of \$0.12/lb proppant and \$1/bbl water (these marginal costs should embody the full marginal cost of the job, and not solely the literal cost of purchasing the raw materials); (d) fixed cost of running the new liner and performing the refrac of \$2M; (e) 10% discount rate; and (f) other reasonable values for taxes, facilities costs, etc. We chose to include the ‘land cost’ in the optimization, and assumed \$15,000/acre.

The ‘objective function’ of the optimization has a big impact on the results. We ran the optimization four times, with different objective functions:

• Maximize DROI, with well spacing set constant at 800 ft (varying only job size)
• Maximize DROI (varying well spacing and job size)
• Maximize NPV/section (varying well spacing and job size)
• Maximize a hybrid objection function weighting both DROI and NPV/section (varying well spacing and job size)

DROI is defined as: (time-discounted revenue minus discounted OPEX) divided by (discounted CAPEX). A DROI of 1.0 breaks even.

The figure below shows the optimal frac design maximizing DROI, enforcing an 800 ft well spacing. The ‘grey’ parts of the fractures show the boundary with the no-flow boundary condition to account for offset wells. The refrac uses 2444 lbs/ft and ~20 bbl/ft.

The cross plot below shows EUR versus total lbs/ft of proppant (800 lbs/ft from the original frac job, plus the lbs/ft from the refrac). EUR reaches the point of diminishing returns at around 3000 lbs/ft.

The crossplot below shows the relationship between DROI and lbs/ft (800 lbs/ft from the original frac job, plus the lbs/ft from the refrac). DROI is maximized at around 3000 lbs/ft, the concentration where the EUR versus job size plot shows the onset of diminishing returns. The maximum DROI is 2.04.

In the next optimization, the code was allowed to vary well spacing in order to maximize DROI. The optimal well spacing is quite wide – 1500-2000 ft. The optimization spaced the wells widely in order to maximize return on investment, even at the cost of reduced NPV and EUR per section. Because of the wide spacing, the job size is huge – 5000 lbs/ft proppant.

This may seem unreasonably large, but if an operator really did want to maximize DROI at the expense of EUR or NPV/section, wide spacing and huge job size might actually be the best design (even though this would be very unconventional). However, in practice, an operator would not want to sacrifice so much on reserves, and this design would not be likely to ever be adopted.

Next, we performed an optimization to maximize NPV/section. The optimal well spacing was 200 ft. The job size on the refrac was 2300 lbs/ft. The NPV/section was \$82M. In contrast, the previous case with 2000 ft spacing had NPV/section of \$28M. DROI was 1.76. The previous case had DROI of 2.14. EUR per 10,000 ft lateral was 310k; the previous case had EUR of 857k per 10,000 ft lateral. However, on a per section basis, the 200 ft spacing design had EUR/section that was 3.6x higher.

200 ft spacing might sound excessively tight, but there are, in fact, some very successful companies using such tight well spacings in the Eagle Ford!

Finally, we ran an optimization with a blended objective function that was 50% DROI and 50% NPV/section. In this case, the optimal well spacing was 400 ft, with 2060 lbs/ft in the refrac. DROI was 1.93, and NPV/section was $64M.

 

The cross-plot of DROI versus NPV/section (from the hybrid objective function optimization) shows that above around $30M/section, DROI decreases as well spacing is tightened to increase NPV/section.

This plot of NPV/section versus well spacing shows how tightening spacing leads to higher NPV/section.

This plot of DROI versus spacing shows that the tightest well spacings (which maximize NPV/section) do not maximize DROI. DROI reaches a maximum at an intermediate spacing, and then decreases beyond that.

The table below summarizes the results from the four optimization runs. The values of well spacing, job size, NPV/section, EUR/section, and DROI vary widely. This exercise demonstrates the tradeoffs available to operators – they can increase EUR and NPV/section by tightening spacing, but at the cost of return on investment. The key for each company is to identify the sweet spot for their own specific priorities.

In an upcoming blog post, I will explore more ideas for the optimization ‘objective function.’ Ultimately, a company’s goal is to maximize its terminal value, given its land position and capital budget. This requires optimization to go beyond the design of a single pad, and to think about the long-term development strategy of the entire asset. Given the major impact of the ‘objective function’ on the design selection, it deserves careful consideration.

 

Recap

This case study demonstrates a refrac modeling and history matching workflow, followed by economic optimization of design. The key takeaways:

• History matching refrac data requires calibration of the degree of ‘crossflow’ between the new perfs and the original fractures. Some amount of fluid pumped back into the preexisting fractures is inevitable and reduces overall efficiency.
• The economic ‘objective function’ has a major impact on the optimization of well spacing and job design. Tighter well spacing trades greater EUR and NPV per section against DROI. Regardless, the job size needs to be dynamically optimized for different values of well spacing.

References

Brinkley, Kourtney, Trevor Ingle, Jackson Haffener, Philip Chapman, Scott Baker, Eric Hart, Kyle Haustveit, and Jon Roberts. 2021. An Eagle Ford Case Study: Monitoring Fracturing Propagation Through Sealed Wellbore Pressure Monitoring. Paper SPE-204140-MS presented at the SPE Hydraulic Fracturing Technology Conference, virtual.

Kang, Charles A., Mark W. McClure, and Somasekhar Reddy, Mariyana Naidenova and Zdravko Tyankov. 2022. Optimizing Shale Economics with an Integrated Hydraulic Fracturing and Reservoir Simulator and a Bayesian Automated History Matching and Optimization Algorithm. Paper SPE-209169-MS presented at the SPE Hydraulic Fracturing Technology Conference, The Woodlands, TX.

Maity, Debotyam, and Jordan Ciezobka. 2020. A Data Analytics Framework for Cored Fracture Imaging and Novel Characterization Workflow – Application on Samples from Hydraulic Fracturing Test Site HFTS in the Midland Basin. Paper SPE-199754-MS presented at the SPE Hydraulic Fracturing Technology Conference and Exhibition, The Woodlands, TX.

McClure, M., C. Kang, S. Medam, C. Hewson, and E. Dontsov. 2021. ResFrac Technical Writeup. https://arxiv.org/abs/1804.02092.

Rezaei, Ali, Giorgio Bornia, Mehdi Rafiee, Mohamed Soliman, and Stephen Morse. 2018. Analysis of refracturing in horizontal wells: Insights from the poroelastic displacement discontinuity method. International Journal for Numerical and Analytical Methods in Geomechanics, 1-22.