The synergy among multiple disciplines, which is achieved based on effective use of available data, is essential to improve the understanding of the physical phenomena behind reservoir responses from current practices of hydraulic fracturing, especially now that the exploitation of unconventional reservoirs is experiencing a challenging time due to oil and gas low-price market conditions. Nevertheless, current market restrictions can be taken as a valuable opportunity to focus on studies that evaluate the status of field pads and design new strategies to maximize returns. This study will present an extension of the integrated project presented by Morales et al.^[1] with an enhanced assessment of the effectiveness of propped and unpropped discrete fracture characterization along complex hydraulic fractures as a result of the synergy of the mentioned integrated characterization project.

The proposed methodology of this research considers applying a non-intrusive EDFM to tune the complex fractures with numerical reservoir simulation. The EDFM method can accurately and efficiently model any complex fracture geometries in reservoir simulation^{[2⇓⇓⇓-6]}. The non-intrusive EDFM method can be employed with any conventional reservoir simulator to model any complex fractured system, by which fracture geometries and attributes can be straightforwardly and accurately embedded in the configured matrix blocks, without the need to migrate to local grid refinement (LGR) or unstructured gridding technique^{[2,7 -8]} (These modeling methods often demand large computational resources for a single well model^[9], which might complicate more the analysis of these very complex fracture systems).

This study is a continuation of Morales et al.^[1], from which we add a novel workflow and resume the field application in Eagle Ford to model and evaluate representative complex hydraulic fractures with the EDFM method. In the first part, we summarized the work of Morales et al.^[1]whose data integration (from seismic interpretation, core and logging data, chemical tracers, fiber optics, pressure-volume-temperature (PVT) analysis, and production analysis) into numerical reservoir modeling. Then, the EDFM method and the new algorithm were used to enhance the research results of Morales et al.^[1] to assess the complex fracture characterization.

The methodology proposed in this study expands and enhances the workflow presented by Morales et al.^[1], regarding a previous field development in Eagle Ford through an integrated characterization project from multiple data sets of different disciplines. Our workflow focuses on the final stage of the characterization process, which refers to modeling the complex hydraulic fractures with EDFM, making the process more efficient in terms of time and handiness since this method has no complication for gridding around fractures, and thereby, large number of complex fractures can be conveniently simulated. The ground upon which Morales et al.^[1] investigation stands is to overcome complex data acquisition synergy after compiling different workflows as seen in Fig. 1. First, microseismic was acquired on three wells to visualize how fracture geometries are disturbed by different completion design. Also, cement bond logs (CBL), ultra-sonic image tool (USIT), and gyros were taken to verify zonal isolation along the wellbores. Multiple facture calibration tests (FCT) were collected in the study area at different target intervals to calibrate vertical stress profiles and examine reservoir pressures. Hydraulic fracture extent and drainage trajectory were defined by water and oil-soluble tracers. Surface pressure gauges at the parent wells aid to characterize hydraulic fracture hits and validate microseismic events (MSE) extent. Furthermore, core and quad-combo were taken in the area to calibrate petrophysical and geomechanical properties. An equation of state (EOS) from oil samples was derived for fluid characterization and reservoir simulation. Lastly, following one year of production, a temporary fiber optic production log was run to determine cluster contribution based on completion design. The independent datasets were integrated on a common commercial software platform for geomodel creation, hydraulic fracture simulations and reservoir simulations (Fig. 1).

In that sense, all the parameters of the integrated model were reviewed in order to validate and extract the geological, petrophysical, PVT, rock physics functions, and fracture intrinsic properties. Next, the complex fracture geometries and properties were extracted from pseudo 3D fracture propagation models and transformed to EDFM input data. In order to do this, the complex fractures were screened by their propped-status and assigned with different fracture conductivity values. The final fracture setting was coupled through EDFM into a commercial reservoir numerical simulator. Moreover, it was relevant to assign different compaction curves for each propped-type fracture group in order to capture the change in the rock compressibility through time due to pressure depletion and stress-dependent mechanical behavior of the permeability along the fractures. Finally, a comprehensive characterization led to achieve satisfactory history match, whose results were employed to estimate ultimate recoveries after 30 years.

The present study is performed in a region of Eagle Ford of south Texas, USA, which is an unconformity-bounded depositional sequence between Cretaceous Buda Formation (bottom) and the Austin Chalk (top). The classic Eagle Ford depositional sequence consists of a retrogradation lower member and a progradational upper member, both of which are separated at or near the maximum flooding surface.^[10-11]In the study area, most of the upper member of the Eagle Ford is truncated beneath the unconformity. As a consequence, core measurements were used to calibrate the petrophysical and geomechanical model. Logging methods used in petrophysical and geomechanical analysis include diameter, gamma ray, resistivity, neutron density, acoustic, and porosity logs^[1]. Parameters in geomechanical analysis include Biot coefficient, elastic modulus, pore pressure and Poisson's ratio^[1]. Heterogeneity is characterized by four electrofacies models (limestone, organic-rich marlstone, organic-poor marlstone, and high-clay content) on the core well after using a combination of normalized gamma ray, deep resistivity, and compressional sonic logs.

This study is applied on a pad with three wells (Well-8, Well-9, and Well-10) in a wine-rack staggered configuration to evaluate complex fracture hits targeting upper and lower Eagle Ford. Lateral segments were drilled at approximately N299° deviation azimuth, which is coherent with the minimum horizontal stress direction in the region. The average horizontal spacing among these wellbores is 100m (330 ft), and vertical spacing reaches up to 18.3 m (60 ft).

Furthermore, the appraised wellbores included multiple unique completion designs established by the operator, as shown in Fig. 2. SW, XL, HY are slick water, crosslink, and hybrid respectively in Fig. 2. Taking HY-8-20- 2000RR as an example, it implies completion design with hybrid fluid, clusters count 8, cluster spacing 6 m (20 ft), and proppant loading per cluster 907 kg (2000 lb), and RR represents increased pumping rate. As a reference, the operator verified the absence of improper zonal isolation along the laterals. More details about the completion schemes and the field operations that were performed on these three wells were discussed by Morales et al.^[1].

As part of the previous integrated study, four total fracture sets were included in the field DFN model based on height variations from outcrop and type (i.e., mixed vs. closed) from image log. In fact, seven different types of fractures are observed in the borehole FMI log analysis: close (153), conductive (1), fault plane (1), open (11), mixed (63), resistive (49), and shear (4). However, mixed, resistive, and conductive fractures were combined into a set called mixed fractures in the final DFN model. Furthermore, the global P₃₂ (area of the fractures per volume of rock) extrapolation also considered the proportionality of the intensity of each fracture set. Fig. 3 exhibits the seismically derived P32 variation for all fractures in the system.More details of the DFN spatial distribution statistics and its construction process are presented by Morales et al.^[1].

A key component on which this study relies on corresponds to the proper characterization of complex fracture model. In that sense, an Unconventional Fracture Model (UFM)^[12-13] was developed in a commercial software to incorporate the geological heterogeneity, geomechanical data, DFN, and microseismic data. The UFM is an extension of a P3D model with the capability to interact with multiple explicit natural fractures and change fracture direction based on a hydraulic fracture propagation / natural fracture crossing criterion^[14]. For instance, a side view of the simulated hydraulic fractures is presented in Fig. 4. Likewise, Fig. 5 shows a map view of the same simulated stage that now includes a slice of the DFN present in the region. The results which show that complex fracture networks are caused by the interaction between hydraulic fractures and DFN distribution^[15], will be later coupled with EDFM automatically thanks to algorithms developed for this study.

One fundamental novelty value of this study corresponds to the capability coupling commercial 3D fracture propagation simulator with EDFM. Thanks to a new preprocessing algorithm, EDFM can read automatically the UFM complex geometries and most important fracture attributes and translate them explicitly into any numerical reservoir simulator without the need of re-griding the original simulation mesh. In synthesis, the algorithm reads the complex fractures information such as the length units, the coordinates of each fracture plane in irregular polygon shapes, the scale, the number of fracture planes, the number of points of each fracture polygon, and the unit vectors to honor the spatial distribution of the complex hydraulic fractures. For instance, the hydraulic fracture properties (permeability, aperture, proppant type, compressibility, etc.) are collected for each fracture segment and transformed into EDFM reservoir input data. As a consequence, a key advancement of this study is the capability to couple UFM and EDFM smoothly and efficiently for posterior reservoir simulation of numerous complex fractures through multiple time-effective realizations.

The complex fracture networks and geomodel properties were translated into an equivalent structured model through EDFM methodology to insert the complex UFM geometries without the need of employing expensive unstructured gridding for numerical simulation. The original grid extent was not disturbed for the sector to maintain the representation of the complex fracture networks. The resulting simulation grid is shown in Fig. 6. The total number of the grid cells in the sector model was close to 50000, a size that presents huge advantages for simulation despite running more 100 000 complex fracture segments.

The information below contains the simulation related properties that are employed in this investigation. For all scenarios, the simulation model size is 1083 m×2547 m × 26 m (3550 × 8350 × 84 ft), and the grid cell size for the model is 15 m × 15 m × 6 m (50 × 50 × 21 ft). The simulated production time lasts 210 days. The initial reservoir pressure is 51 MPa (7344 psi), reservoir temperature is 107 °C (225°F), and the reservoir depth ranges from 3011 to 3159 m (9877-10365 ft). The average porosity and permeability of matrix are 5% and 130×10^-9 μm², respectively. We collected the information regarding initial water saturation from the geostatic model, and its average value is found at 6%. Furthermore, the detailed information of the fluid model process is covered in depth in a previous work, and it is not subject of analysis in this study.

For this study, the complex fracture networks were discretized and screened by their propped status into three additional regions with respect to the matrix, they are unpropped, propped and fully propped as showed in Fig. 7. These fracture regions are directly related to their sediment size (propped = 150 μm (100 mesh), and “fully propped” = 550/270 μm (30/50 mesh). In fact, the criterion of determining which fracture planes corresponds to each group is based on the results of the UFM simulation, which considers multiple factors when computing quantitative transport physics such as settling velocity for each proppant loading per cluster, and heterogeneous proppant placement in variable-viscosity fluid systems even though some transport properties depend on the overall proppant volume fraction. Therefore, clustering the fractures by these criteria helped to quantify their effectiveness.

It is well known that fracture conductivities are stress- dependent or pressure-dependent critical variables in the characterization of hydraulic fracture settings^[4,16]. In that sense, and as mentioned earlier in this paper, compaction tables were employed to capture the stress-dependent mechanical behavior of the permeability along the fractures due to pore pressure depletion. These compaction tables compute multipliers on porosity, permeability, and/or transmissibility at specific pressure points on the basis of discrete values lookup, and linear interpolation methods. The fundamental purpose of the use of these rock compaction curves is to define the dynamic compressibility path and fracture conductivity behavior as non-linear relationships between porosity, permeability and/or transmissibility at deterministic pressure values.

For this study, as part of the different modeling tuning parameters prior history matching, an experiment was proposed to evaluate the influence of the rock compaction curves on the simulation response of the studied model. The experiment consisted in contrasting the BHP response after applying three gradient levels of the rock compaction curves: with severe, intermediate rock compaction tables, and one without any rock compaction modification. This compaction severity is defined as how low the permeability multiplier is reduced in respect to the pressure. Also, the experiment is justified due to the need to test alternatives to constrain the pressure response of the model to match lower BHP field-reported values. The results are displayed in Fig. 8.

As expected, the influence of the compaction phenomena and its degrees of severity as a result of fracture conductivity alteration for each wellbore can be clearly observed, specially for wells 8 and 9. The more severe the compaction curve is (lower permeability multiplier in respect to the pressure), the more depleted the BHP will be. The case of Well-10 is particular since in all three compaction cases, the pressure depletion difference is narrower among the cases due to higher production volume with respect to the other wells as seen in the following section.

For this research, the history matching process involved calibrating a series of parameters, including, but not limited to, fracture conductivities, fracture water saturation, and rock compaction curves in different regions of proppant-status intensity. The tunning process of these parameters was enhanced and simplified by the use of the non-intrusive EDFM methodology and the previously mentioned algorithms that allowed to interact and test multiple realizations and sensitivities efficiently among the different stages of the workflow. For this validation process, liquid production rates were input constrain, and BHP was used as the responsive history-matching parameter.

For instance, the massive number of connections among the wellbores and the complex fractured system from the UFM final realization (Figs. 6 and 7) posed the greatest challenge in this dynamic characterization so that the pressure drawdown can be represented accurately. Also, the constant fracture heights from the UFM high-resolution complex fractures that crossed all the layers of the reservoir contributed to this unrealistic bottomhole pressure build-up in the numerical domain. As a solution to this issue, a calibration of the Net to Gross (NTG) value and the fracture conductivities were implemented in order to honor the fracture geometries but to achieve the field-reported BHP depletion. As a result, it can be concluded that the fundamental factor with greater impact in the history matching, was the characterization of the dynamic behavior of the fracture conductivities, which were addressed by the pressure-dependent compaction curves for each propped fractured group as explained earlier in this paper. Detailed matching results are presented in Fig. 9. As a consequence, it is clearly observed that the flowing bottomhole pressure of all three wells achieved a satisfactory match. For instance, an educated inference from eyeball estimation of Fig. 9 is that the wells with better history match are Well-9 and Well-10 because of their more stable BHP response that fits within close ranges from the observed field data.

Furthermore, one of the novelties of this study corresponds to the capability of the developed algorithms to generate 3D visualizations of pressure distribution along the multiple planes of the complex fractures thanks to the non-intrusive EDFM (Fig. 10). Furthermore, the proposed workflow includes the capability to isolate the visualization of the fractures from the rest of the field or the matrix blocks. This study reveals that even though the fully-propped fractures region is the first zone to deplete (Fig. 10b, 10d), the other propped fractures of region 2 can contribute to the effective drainage of the complex fractures at later time (Fig. 10f).

As expected, more complex fractures are connected in the fully propped regions, generating more well interference and higher pressure drop. Hence, pressure is lower through time in the fractures closer to the wellbore or in the fully-propped ones because these fractures have higher dimensionless conductivity. This higher dimensionless fracture conductivity indicates excessive fracture efficiency and proposes that the fracture capability to deliver reservoir fluids to the wellbore is more efficient than the formation's capability to deliver the same fluids to the fracture. Nevertheless, we can also identify unpropped regions along some parts of the wellbores. One of the possible reasons behind unpropped regions along the well is that those regions correspond to completions schemes with tighter cluster spacing, which can cause more fractures to grow up and down through stress barriers^[17], resulting in more proppant being placed out of the target zone.

As discussed by Fowler et al.^[18] permeability estimation is a key parameter for the optimization of well completions designs. In that sense, this study considers flow to be strongly influenced by the distribution of proppant and conductivity in the fracture. Moreover, the fracture characterization in this study contemplates that conductivity evolves over time as a function of rock compaction^[19].

Modeling subsurface scenarios with complex fractures is essential to anticipate field practices outcomes^[20]. Also, the engineering human talent should be complemented by effective technology. Therefore, newer EDFM applications can guide operators on how to estimate reliable estimated ultimate recovery (EUR), tailor preventive treatments, predict fracture hits, and redesign child wells forecasting prospective interactions. The modeling results are consistent with the data synergy after careful evaluation.

The performance of EDFM method seems quite convenient in terms of computation time for translating field-scale UFMs with P3D propagation of complex fractures into reservoir numerical simulation. The processes of high-resolution complex-fracture modeling in a reservoir-scale simulation are often extremely time-consuming when working with numerous fracture segments in unstructured gridding. In that sense, significant reduction in the computation time helped the ability to perform more sensitivity studies to calibrate the model. Using the EDFM method, a single case of three wells was able to run in less than an hour, while unstructured gridding cases took more than 5 hours in a standard eight-processor machine.

The pressure depletion can be observed on only propped and fully-propped fractures in the model of this study that employs EDFM method. The propped regions exhibit larger effective hydraulic fracture surface area and more connections to the complex natural fracture network. Finally, the pressure depletion and drainage area results of this study aided to qualify the best completion design when compared to the multiple completion schemes that were used in these three wellbores by the operator. In that sense, two types of fluid with higher oil recovery excelled, and they are the crosslink (XL) in 8 clusters separated 9 m (30 ft) with 544 kg (1200 lb.)/cluster of proppant loading, and the hybrid fluid in 8 clusters separated 6 m (20 ft), and with 907 kg (2000 lb.)/cluster of proppant.

After consistent history match was achieved, production forecast is performed in order to estimate ultimate recoveries over 5, 10, and 30 years. Pressure was constrained to 3.4 MPa (500 psi), and liquid flow rates were forecasted. Detailed results are presented in Fig. 11. It can be inferred from the forecast plots and 3D visualizations that Well-8 is underperforming when compared to the other two wells which share a similar recovery. This underperformance is directly related to the reduced effectiveness of the completions closer to the heel, which triggered unpropped fractures regions along the wellbore in those zones. As a result, this condition reduced the well drainage area, leaving hydrocarbons trapped in the shale for this well.

The novel feature of this study is related to optimizing the complex-fracturing characterization through EDFM in comprehensive unconventional reservoir models by coupling a calibrated geological model, with fracture diagnostic tools, UFM simulation, and production data via faster and friendlier algorithms for further evaluations of multiple field practices.

Optimizing the simulation time of very extensive complex fractures could help to explore different hydraulic communication scenarios instead of using high-computing demanding frameworks. Generating a robust geological model and calibrating the hydraulic fracture propagation model based on the heterogeneous fracture designs along the wellbores trigger the need to validate how influential is the proppant intensity on effective surface area as presented in the results of this study. Effectiveness of modeling complex fractures was constrained by propped high/poor intensity which led to larger/smaller fracture surface area in the network. A fundamental parameter used to achieve satisfactory history match is the non-linear dynamic characterization of pressure-dependent fracture conductivities along the different regions of the fractured system.

We thank Chesapeake Energy for permission to publish this work, and to acknowledge Chesapeake Energy and Sim Tech LLC for providing the data and the software for this investigation. A special thanks to Adrian Morales who contributed technical insights and shared all the information from his previous work.