This study introduces a methodology for estimating uncertainty in production of new shale wells. The methodology combines geostatistical modeling and machine learning and accounts for geological uncertainty. Our approach improves uncertainty quantification by merging local and global trends.
A functional random forest regression model is trained, connecting completion parameters and geological parameters with production profiles. The geological parameters are simulated on a two-dimensional areal grid using Sequential Gaussian Simulation with information from pilot wells. Production profiles are generated for each cell in the simulation grid based on optimal completion and simulated geological parameters. These machine learning-based realizations account for the global trends. We account for local effects by using cokriging to merge total production from nearby wells with the production from machine learning-based realizations.
We test the methodology on a dataset from the Eagle Ford formation. The result of the study is a map of the play highlighting the most probable production (P50) for different areas and associated risk (P90–P10). The resulting map allows us to rank locations for new wells for drilling. The proposed methodology provides a first estimate, and a more detailed data investigation is required to sanction a new well in a particular location.
Jew, Adam (SLAC National Accelerator Laboratory) | Li, Qingyun (Stanford University) | Cercone, David (National Energy Technology Laboratory) | Brown, Gordon E. (Stanford University) | Bargar, John (SLAC National Accelerator Laboratory)
Barite scale formation in unconventional systems can severely reduce production rates and recovery. Barite in drilling mud (DM) is highly susceptible to dissolution (up to 20% of the total barite) when it interacts with a 15% hydrochloric acid spearhead leading to significant barite scaling as solution pH rises to circumneutral levels. In order to maintain the low pH spearhead that is favored by most operators, a new acid spearhead formulation was developed to stabilize the barite in place, closer to the borehole where open volume is large. This approach is dramatically different than typical scale mitigation practices in which precipitation from solution is targeted for control; but with unsatisfactory results. By using sulfuric acid instead of hydrochloric acid (equivalent acidity), the amount of Ba and SO4 leached from DM was reduced by > 850-fold. This major reduction in barite dissolution translates into a significant reduction in the dispersal of these mineral-forming aqueous species in the hydraulically fractured shale matrix. As a result, reprecipitation of barite scale in the fracture network is minimized. Because some unconventional systems have high concentrations of calcite that dissolves upon interaction with the acid spearhead, gypsum can precipitate due to excess Ca2+ derived from the calcite and its interaction with the additional SO42− introduced from the new acid spearhead formulation. To mitigate this potential scale formation, Na-citrate was added to chelate Ca2+ and inhibit gypsum formation. When only H2SO4 was used, Na-citrate was not able to completely inhibit gypsum formation in high carbonate systems. However, by adjusting the H2SO4/HCl ratios and including Na-citrate, both barite and gypsum scale in the stimulated rock volume can be mitigated. These findings suggest that an incremental new best-practice is to stabilize the Ba/SO4 in place rather than attempt to control precipitation from solution.
Barite is one of the most important and problematic types of mineral scale in unconventional shale systems. Depending on the mineralogy of the deposit and solution pH, the amount and location of barite scale can differ (Figure 1). In systems where the generation of secondary porosity is minimal because of low concentrations of carbonate minerals, the precipitation of barite can reduce permeability by > 15% (Alalli et al., 2018). Efforts to mitigate barite precipitation focus primarily on either keeping barite nucleation from occurring or altering the surface properties of newly formed barite crystals to limit adhesion to shale and pipe surfaces (Antonious, 1996; Jones et al., 2003; Shen et al., 2009; Mavredaki et al., 2011; Yan et al., 2015; Bhandari et al., 2016; He and Vidic, 2016; Yan et al., 2017). These methods have limited success due to the intentional introduction of high quantities of barite in unconventional systems through drilling mud (DM). Previous work has shown that dissolution of barite in DM due to the injection of chemicals into the subsurface, primarily the 15% HCl acid spearhead, releases both Ba2+ and SO42− into the system that can precipitate as barite in fractures and shale matrices (Jew et al., 2018). On a per-mass basis, the barite released from the DM by the acid spearhead is significantly higher than that from the shale itself, indicating that the DM is most likely the dominant Ba/SO4 source in these systems. Because of this underappreciated Ba/SO4 source in unconventional systems, we developed a new model for barite in the subsurface (Figure 2). Because Ba and SO4 are introduced into the system before stimulation fluid injection, any barite dissolved from the DM by the acid spearhead will be transported into stimulated rock volume (SRV) by the slickwater/slurry where precipitation will occur, necessitating a new approach to scale mitigation.
Li, Qingyun (SLAC National Accelerator Laboratory / Stanford University) | Jew, Adam (SLAC National Accelerator Laboratory) | Cercone, David (National Energy Technology Laboratory) | Bargar, John (SLAC National Accelerator Laboratory) | Brown, Gordon E. (SLAC National Accelerator Laboratory / Stanford University) | Maher, Katherine (Stanford University)
Laboratory experiments have shown that hydraulic fracturing fluids (HFF) can chemically interact with iron(Fe)-bearing minerals in shale, releasing Fe(II) which is then oxidized to form Fe(III)-(hydr)oxide scale. The Fe(III)-(hydr)oxide scale can occlude pore space and reduce oil and gas production in wells. Our previous experimental studies show that Fe(III)-(hydr)oxides can precipitate even under acidic conditions where Fe(II) oxidation is unexpected. This is due to bitumen that is extracted from shale by chemical additives in HFF, which can aid in Fe(II) oxidation and lead to the precipitation of Fe(III)-(hydr)oxides. In this numerical modeling study, we built two geochemical models to simulate our experimental observations. One model was developed to construct the rate law of Fe(II) oxidation in the presence of bitumen in a shale-free system, while the other model was used to understand the chemical reaction network and quantify the impact of bitumen on iron oxidation/precipitation during shale-HFF interactions. Our modeling results have shown that in both high- and low-carbonate shale systems, the presence of bitumen can increase the rate of iron oxidation/precipitation by more than an order of magnitude compared to bitumen-free systems. In addition, the availability of dissolved oxygen to pyrite grains is critical for Fe dynamics during shale-HFF interactions. The chemical reaction network obtained from this study, along with the bitumen-aided Fe(II) oxidation rate law, pave the way for future modeling studies on chemical reactions during hydraulic fracturing and their influence on formation damage and lost hydrocarbon production.
Scale formation during unconventional stimulation can reduce porosity and permeability which, in turn, can affect hydrocarbon flow in oil and gas wells. Operators add iron (Fe) control agents such as citrate and methanol to the hydraulic fracturing fluid (HFF) to prevent the deposition of Fe(III)-(hydr)oxide (Fe(OH)3) minerals in the pore space of shale (FracFocus). Our previous studies have shown that even with Fe control agents, Fe(OH)3 can still precipitate (Harrison et al. 2017, Jew et al. 2017). The Fe in these Fe(III)-bearing precipitates comes primarily from dissolution of pyrite (FeS2), a common mineral in shale, during shale-HFF interactions. In unconventional reservoirs, Fe is released from corrosion of downhole steel in addition to dissolution of native pyrite due to injection of a 15% hydrochloric acid spearhead, resulting in more severe scaling. During pyrite dissolution and pipe corrosion, Fe is released to fluid in the Fe(II) oxidation state, and is further oxidized by dissolved oxygen (O2) in the fluid to Fe(III), which readily precipitates upon acid neutralization. Our previous studies also found that Fe(II) oxidation can be aided by aqueous bitumen leached from shale by organic compounds in the fracturing fluid (Jew et al. 2017). The effect of bitumen on Fe(II) oxidation is enormous, especially under acidic conditions, where bitumen was found to promote Fe(II) oxidation and Fe(OH)3 precipitation that otherwise would not occur. This phenomenon is important, because in most hydraulic fracturing operations, the acid spearhead is the first chemical injected down the borehole for cleaning purposes. Subsequent injections of slickwater/slurry allow for the acid to be pushed into the stimulated rock volume (SRV). Under the acidic conditions along the fluid flow pathways, the formation of Fe(III)-(hydr)oxide scale depends on not only the dissolved O2 concentration in solution but also bitumen-aided Fe(II) oxidation.
In this study we explore the use of multilevel derivative-free optimization for history matching, with model properties described using PCA-based parameterization techniques. The parameterizations applied in this work are optimization-based PCA (O-PCA) and convolutional neural network-based PCA (CNN-PCA). The latter, which derives from recent developments in deep learning, is able to represent accurately models characterized by multipoint spatial statistics. Mesh adaptive direct search (MADS), a pattern search method that parallelizes naturally, is applied for the optimizations required to generate posterior (history matched) models. The use of PCA-based parameterization reduces considerably the number of variables that must be determined during history matching (since the dimension of the parameterization is much smaller than the number of grid blocks in the model), but the optimization problem can still be computationally demanding. The multilevel strategy introduced here addresses this issue by reducing the number of simulations that must be performed at each MADS iteration. Specifically, the PCA coefficients (which are the optimization variables after parameterization) are determined in groups, at multiple levels, rather than all at once. Numerical results are presented for 2D cases, involving channelized systems (with binary and bimodal permeability distributions) and a deltaic-fan system, using O-PCA and CNN-PCA parameterizations. O-PCA is effective when sufficient conditioning (hard) data are available, but it can lead to geomodels that are inconsistent with the training image when these data are scarce or nonexistent. CNN-PCA, by contrast, can provide accurate geomodels that contain realistic features even in the absence of hard data. History matching results demonstrate that substantial uncertainty reduction is achieved in all cases considered, and that the multilevel strategy is effective in reducing the number of simulations required. It is important to note that the parameterizations discussed here can be used with a wide range of history matching procedures (including ensemble methods), and that other derivative-free optimization methods can be readily applied within the multilevel framework.
Accurate and robust well modeling is essential for performing reservoir simulations of practical interest. The Multi-Segment well (MSWell) model is able to describe the well topology and accurately represent the multiphase multicomponent flow and transport behavior in the wellbore. The fully coupled method (FC) has been developed and widely applied on coupled reservoir and MSWell modeling due to its unconditional stability and consistent implementation. A local well solver can be applied to provide a better nonlinear precondition for MSWell variables in order to accelerate the nonlinear convergence of the FC method.
However, solving the coupled MSWell and reservoir model in a fully implicit scheme can still present limitations on some practical applications. First, the well or surface facility solver can be separate from the existing reservoir simulator, making it challenging to employ the fully implicit method. Second, complex linear and nonlinear solvers need to be designed to pair the specific wells and reservoir models. These solvers have to account for the different flow characteristics and discretization domains between reservoir and MSWell. A sequential coupling scheme can become preferable in such situations.
Sequential fully Implicit method (SFI) splits the fully coupled reservoir and MSWell equations into two parts and solves them sequentially. In spite of accomplishing an implicit coupling in a sequential scheme, SFI suffers the slow outer loop convergence rate especially when reservoir is strongly coupled with the wells, which is very often the case. The slow convergence is caused by the linear convergence rate of the fix point iteration used in the SFI. Here, we developed a sequential implicit Newton's method (SIN) for coupled MSWells. SIN incorporates a Newton update at the end of each sequential step to achieve a quadratic convergence of outer iterations, while require a limited extra computational cost. Numerical results show that SIN attains comparable nonlinear Newton iterations with the FC in the coupled heterogeneous reservoir and complex MSWell problems.
Pan, Huanquan (Stanford University) | Imai, Motonao (Japan Oil, Gas and Metals National Corporation/Waseda University) | Connolly, Michael (Stanford University) | Tchelepi, Hamdi (Stanford University)
Robust and efficient multiphase flash calculations are crucial in compositional and thermal simulations for complex fluid systems in which three or four phase may co-exist. Solution of the Richford-Rice (RR) equations is an important operation in the multiphase flash. The Newton method generally does not converge during solution of the RR equations unless very good initial values are provided. In this paper, the solution of the RR equations is formulated as a minimization of a convex function problem. For the first time, we use a trust-region (TR) method to solve the RR equations through minimization of the convex function. The Hessian matrix of the convex function is always positive-definite, and the TR-based solver guarantees convergence. The key to successful implementation is to determine the relaxation parameter in the Newton update. We select this relaxation parameter to meet the boundary of the objective function and to ensure an adequate step length. We tested the RR solver for three and four phase RR problems in the construction of phase diagrams. The test cases are representative of complex fluid systems encountered in enhanced oil recovery, including injection of CO2 into low temperature reservoirs and steam injection into heavy oil reservoirs at elevated temperatures. We performed tens of millions of multiphase flash computations, the results of which reveal our RR solver to be robust, efficient, insensitive to initial values, and capable of handling negative phase amounts. We also evaluated the effect of the initial values on convergence and recommend methods to estimate the initial values in our RR solver. In summary, our RR solver greatly improves the multiphase flash calculations and strengthens the coupling of phase equilibrium calculations to the governing equations in multiphase compositional and thermal simulation.
A reduced-order modeling framework is developed and applied to simulate coupled flow-geomechanics problems. The reduced-order model is constructed using POD-TPWL, in which proper orthogonal decomposition (POD), which enables representation of the solution unknowns in a low-dimensional subspace, is combined with tra jectory piecewise linearization (TPWL), where solutions with new sets of well controls are represented via linearization around previously simulated (training) solutions. The over-determined system of equations is pro jected into the lowdimensional subspace using a least-squares Petrov-Galerkin procedure, which has been shown to maintain numerical stability in POD-TPWL models. The states and derivative matrices required by POD-TPWL, generated by an extended version of Stanford's Automatic-Differentiation-based General Purpose Research Simulator, are provided in an offline (pre-processing or training) step. Offline computational requirements correspond to the equivalent of 5-8 full-order simulations, depending on the number of training runs used. Runtime (online) speedups of O(100) or more are typically achieved for new POD-TPWL test-case simulations. The POD-TPWL model is tested extensively for a 2D coupled problem involving oil-water flow and geomechanics. It is shown that POD-TPWL provides predictions of reasonable accuracy, relative to full-order simulations, for well-rate quantities, global pressure and saturation fields, global maximum and minimum principal stress fields, and the Mohr-Coulomb rock failure criterion, for the cases considered. A systematic study of POD-TPWL error is conducted using various training procedures for different levels of perturbation between test and training cases. The use of randomness in the well bottom-hole pressure profiles used in training is shown to be beneficial in terms of POD-TPWL solution accuracy. The procedure is also successfully applied to a prototype 3D example case.
Understanding the effect of injected water salinity is becoming crucial, as it has been shown to have a strong impact on the efficiency of oil recovery process. Various experiments have concluded that carbonate wettability is altered when the water ionic-composition is changed. In this work, a numerical investigation of an oil blob mobilized by water is conducted inside a single pore. The presented model studies the synergy effect of multiphase flow and water salinity at the pore level.
To model the multiphase flow at the pore-scale, the full hydrodynamic Navier-Stokes equations are solved using direct numerical simulation (DNS). The effect of brine ionic-composition is examined through the electric double layer effect. Experimental zeta potential values, published in the literature, of crude oil/water and water/carbonate interfaces have been employed in the model, which capture the water salinity effect.
The simulation results show that the water wetting film surrounding the oil-droplet collapses to an adsorbed nanometer water layer when high salinity water is used. As a result, a large pressure gradient is required to mobilize the oil inside the pore and overcome the attractive surface forces between the oil/water and water/carbonate interfaces. For low-salinity injected water, the carbonate surface becomes more water-wet. The wetting film surrounding the oil blob becomes stable due to the repulsive electric double layer force. Therefore, less energy is required to mobilize the oil blob inside the pore compared to high water salinity. The effect of solid roughness and injected water flow rate are also studied, which show to have a strong impact on the oil displacement efficiency.
The novelty of the numerical method lies in efficiently capturing the nanoscale effect of the electric double layer in pore-scale multiphase flow at the microscale. The simulation results provide fundamental insights on the efficiency of low-salinity waterflooding at the pore level.
The critical gas saturation in permeable sands was studied as a function of depletion rate and the presence of an aqueous phase as the major experimental variables. Voidage-replacement ratios (VRR = injected volume/produced volume) less than 1 were used to obtain pressure depletion with active water injection. Three different live crude oils were considered. Two of the oils are viscous Alaskan crudes with dead-oil viscosities of 87.7 and 600 cp, whereas the third is a light crude oil with a dead-oil viscosity of 9.1 cp. The critical gas saturation for all tests ranged from 4 to 16%. These values for critical gas saturation are consistent with the finding that the gas phase displayed characteristics similar to those of a foamy oil. For a given oil and depletion rate, the critical gas saturation was somewhat larger for VRR = 0 than it was for VRR = 0.7. The oil recovery correlates with the critical gas saturation (i.e., for a given VRR, tests exhibit greater oil recovery when the critical gas saturation is elevated). For the conditions tested, there was not a strong correlation of critical gas saturation over more than two orders of magnitude of the rate of pressure depletion, for a given VRR. Such behavior might be consistent with theoretical studies reported elsewhere that suggest that the critical gas saturation is independent of the pressure-depletion rate when the rate of depletion is small.
This paper presents an analysis of the interactions between stimulation design and two important geomechanical effects: the variation of least principal stress (Shmin) between lithological layers and the stress shadow effect that arises from simultaneously propagating adjacent hydraulic fractures. To demonstrate these interactions, hydraulic fracture propagation is modeled with a 5-layer geomechanical model representing an actual case study. The model consists of a profile of Shmin measurements made within, below and above the producing interval. The stress variations between layers leads to an overall upward fracture propagation and proppant largely above the producing interval. This is due to interactions between the pressure distribution within the fracture and the stress contrast in the multiple layers. A sensitivity study is done to investigate the complex 3-D couplings between geomechanical constraints and well completion design parameters such as landing zone, cluster spacing, perforation diameter, flow rate and proppant concentration. The simulation results demonstrate the importance of a well characterized stress stratigraphy for prediction of hydraulic fracture characteristics and optimization of operational parameters.