Summary In this study, we demonstrate geomechanical modeling with fully automatic parameter calibration to estimate the full geomechanical stress fields of a prospective US CO2 storage site, based on sparse measurement data. The goal is to compute full stress tensor field estimates (principal stresses and orientations) that are maximally compatible with observations within the constraints of the model assumptions, thereby extending point-wise, incomplete partial stress measurement to a simulated full formation stress field, as well as a rough assessment of the associated error. We use the Perch site, located in Otsego Country, Michigan, as our case study. Input data consists of partial stress tensor information inferred from in-situ borehole tests, geophysical well logs and processing of seismic data. A static earth model of the site was developed, and geomechanical simulation functionality of the open-source MATLAB Reservoir Simulation Toolbox (MRST) used to model the stress field. Adjoint-based nonlinear optimization was used to adjust boundary conditions and material properties to calibrate simulated results to observations. Results were interpreted through a Bayesian framework. The focus of this article is to demonstrate how the fully automatic calibration procedure works and discuss the results obtained but does not attempt a detailed analysis of the stress field in the context of the proposed CO2 storage initiatives. Our work is part of a larger effort to non-invasively determine in-situ stresses in deep formations considered for CO2 storage. Guided by previously published research on geomechanical model calibration, our work presents a novel calibration approach supporting a potentially large number of linear or nonlinear calibration parameters, in order to produce results optimally agreeing with available measurements and thus extend partial point-wise estimates to full tensor fields compatible with the physics of the site.

Abstract Radio frequency (RF) heating is recognized as a technique having the potential to thermally enhance remediation of hydrocarbon-impacted soil. RF heating delivers electromagnetic (EM) power to a targeted body of soil, resulting in an increased soil temperature that enhances the in-situ remediation processes such as biodegradation. Antennas are placed either on the ground or installed in the soil near the ground surface. The antennas operate in the hundreds of kHz to MHz range. To model the RF heating process, we successfully coupled a reservoir simulator with a 3-dimensional (3D) EM solver to evaluate the ability of RF technology to heat soil in situ. The coupled reservoir/EM simulator solves the EM fields and associated heating for a heterogeneous reservoir or soil volume in the presence of multiple antennas. The coupling was accomplished through a flexible interface in the reservoir simulator that allows the runtime loading of third-party software libraries with additional physics. This coupled workflow had been previously used for studying RF heating for heavy oil recovery (Li 2019). An RF heating simulation case study was performed in support of a soil remediation field test designed to demonstrate the ability to heat soils using EM energy. The study included field test data analysis, simulation model building, and history matching the model to test data. Results indicate, on average, the soil was heated ∼2-3°C above the initial formation temperature after approximately two days (52 hours) of RF heating. We found that the RF heating was local, and our simulation model, after tuning input parameters, was able to predict a temperature profile consistent with the field test observations. With properly designed RF heating field pilots and tuning of EM and reservoir parameters in simulation models, the coupled reservoir/EM simulator is a powerful tool for the calibration, evaluation, and optimization of RF heating operations.

Abstract The process of extracting coalbed methane (CBM) is not only of significance for unconventional energy supply but also important in mine safety. The recent advance in fracking techniques, such as carbon dioxide (CO2) fracking, intensifies the complexity of stimulated coalbeds. This work focuses on developing a fully coupled multidomain model to describe and get insight into the process of CBM extraction, particularly from those compound-fractured coalbeds. A group of partial differential equations (PDEs) are derived to characterize gas transport from matrix to fractures and borehole. A stimulated coalbed is defined as an assembly of three interacting porous media: matrix, continuous fractures (CF) and radial primary hydraulic fracture (RF). Matrix and CF constitute a dual-porosity-dual-permeability system, while RF is simplified as an 1-D cracked medium. These media further form three distinct domains: non-stimulated reservoir domain (NSRD), stimulated reservoir domain (SRD) and RF. The effects of coal deformation, heat transfer, and non-thermal sorption are coupled into the model to reflect the multiple processes in CBM extraction. The finite element method is employed to numerically solve the PDEs. The proposed model is verified by comparing its simulation results to a set of well production data from Southern Qinshui Basin in Shanxi Province, China. Great consistency is observed, showing the satisfactory accuracy of the model for CBM extraction. After that, the difference between various stimulation patterns is presented by simulating the CBM extraction process with different stimulation patterns including (1) unstimulated coalbed; (2) double-wing fracture + NSRD; (3) multiple RFs + NSRD; (4) SRD + NSRD and (5) multiple RFs + SRD + NSRD. The results suggest that Pattern (5) (often formed by CO2 fracking) boosts the efficiency of CBM extraction because it generates a complex fracture network at various scales by both increasing the number of radial fractures and activating the micro-fractures in coal blocks. Sensitivity analysis is also performed to understand the influences of key factors on gas extraction from a stimulated coalbed with multiple domains. It is found that the distinct properties of different domains originate various evolutions, which in turn influences the CBM production. Ignoring thermal effects in CBM extraction will either overestimate or underestimate the production, which is the net effect of thermal strain and non-isothermal sorption. The proposed model provides a useful approach to accurately evaluate CBM extraction by taking the complex evolutions of coalbed properties and the interactions between different components and domains into account. The importance of multidomain and thermal effects for CBM reservoir simulation is also highlighted.

Abstract Acid fracturing technique is widely applied to stimulate the productivity of carbonate reservoirs. The acid-fracture conductivity is created by non-uniform acid etching on fracture surfaces. Heterogeneous mineral distribution of carbonate reservoirs can lead to non-uniform acid etching during acid fracturing treatments. In addition, the non-uniform acid etching can be enhanced by the viscous fingering mechanism. For low-perm carbonate reservoirs, by multi-stage alternating injection of a low-viscosity acid and a high-viscosity polymer pad fluid during acid fracturing, the acid tends to form viscous fingers and etch fracture surfaces non-uniformly. To accurately predict the acid-fracture conductivity, this paper developed a 3D acid fracturing model to compute the rough acid fracture geometry induced by multi-stage alternating injection of pad and acid fluids. Based on the developed numerical simulator, we investigated the effects of viscous fingering, perforation design and stage period on the acid etching process. Compared with single-stage acid injection, multi-stage alternating injection of pad and acid fluids leads to narrower and longer acid-etched channels.

Abstract Accounting for poro-mechanical effects in full-field reservoir simulation studies and uncertainty quantification workflows is still limited, mainly because of their high computational cost. We introduce a new approach that couples hydrodynamics and poro-mechanics with dual-porosity flow diagnostics to analyse how poro-mechanics could impact reservoir dynamics in naturally fractured reservoirs without significantly increasing computational overhead. Our new poro-mechanical informed dual-porosity flow diagnostics account for steady-state and singlephase flow conditions in the fractured medium while the fracture-matrix fluid exchange is approximated using a physics-based transfer rate constant which models two-phase flow using an analytical solution for spontaneous imbibition or gravity drainage. The deformation of the system is described by the dualporosity poro-elastic theory, which is based on mixture theory and micromechanics to compute the effective stresses and strains of the rock matrix and fractures. The solutions to the fluid flow and rock deformation equations are coupled sequentially. The governing equations for fluid flow are discretised using a finite volume method with two-point flux-approximation while the governing equations for poro- mechanics are discretised using the virtual element method. The solution of the coupled system considers stress-dependent permeabilities for fractures and matrix. Our framework is implemented in the open source MATLAB Reservoir Simulation Toolbox (MRST). We present a case study using a fractured carbonate reservoir analogue to illustrate the integration of poro-mechanics within the dual-porosity flow diagnostics framework. The extended flow diagnostics calculations enable us to quickly screen how the dynamics in fractured reservoirs (e.g. reservoir connectivity, sweep efficiency, fracture-matrix transfer rates) are affected by the complex interactions between poro-mechanics and fluid flow where changes in pore pressure and effective stress modify petrophysical properties and hence impact reservoir dynamics. Due to the steady-state nature of the calculations and the effective coupling strategy, these calculations do not incur significant computational overheads. They hence provide an efficient complement to traditional reservoir simulation and uncertainty quantification workflows as they enable us to assess a broader range of reservoir uncertainties (e.g. geological, petrophysical and hydro-mechanical uncertainties). The capability of studying a much broader range of uncertainties allows the comparison and ranking from a large ensemble of reservoir models and select individual candidates for more detailed full-physics reservoir simulation studies without compromising on assessing the range of uncertainties inherent to fractured reservoirs.

Abstract Complex well model has proved to be important for capturing the full physics in wellbore, including pressure losses, multiphase effects, and advanced device modelling. Numerical instability may be observed especially when the well is produced at a low rate from a highly productive multi-phase zone. In this paper, a new multi-level nonlinear solver is presented in a state-of-the-art parallel complex wellbore model for addressing some difficult numerical convergence problems. A sequential two-level nonlinear solver is implemented, where the inner solver is used to address the convergence in the constraint rate equation, and then the entire complex network is solved using an outer solver. Finally, the wellbore model is coupled with the grid solution explicitly, sequentially, or implicitly. This novel formulation is robust enough to greatly improve the numerical stability due to the lagging in the computation of mixture density in wellbore constraint rate equation and the variation in the fluid composition over Newton iterations in network nonlinear solver. The numerical challenge in the complex well model and the improvement of performance with the new nonlinear solver are demonstrated using reservoir simulation. Models with complex wells running into convergence problems are constructed and simulated. With this novel nonlinear solver, simulation gives much more reliable results on well productions without numerical oscillations and computational cost is much less.

Abstract Alternative to CPU computing architectures, such as GPU, continue to evolve increasing the gap in peak memory bandwidth achievable on a conventional workstation or laptop. Such architectures are attractive for reservoir simulation, which performance is generally bounded by system memory bandwidth. However, to harvest the benefit of a new architecture, the source code has to be inevitably rewritten, sometimes almost completely. One of the biggest challenges here is to refactor the Jacobian assembly which typically involves large volumes of code and complex data processing. We demonstrate an effective and general way to simplify the linearization stage extracting complex physics-related computations from the main simulation loop and leaving only an algebraic multi-linear interpolation kernel instead. In this work, we provide the detailed description of simulation performance benefits from execution of the entire nonlinear loop on the GPU platform. We evaluate the computational performance of Delft Advanced Research Terra Simulator (DARTS) for various subsurface applications of practical interest on both CPU and GPU platforms, comparing particular workflow phases including Jacobian assembly and linear system solution with both stages of the Constraint Pressure Residual preconditioner.

Abstract In many applications the accurate representation of the computational domain is a key factor to obtain reliable and effective numerical solutions. Curved interfaces, which might be internal, related to physical data, or portions of the physical boundary, are often met in real applications. However, they are often approximated leading to a geometrical error that might become dominant and deteriorate the quality of the results. Underground problems often involve the motion of fluids where the fundamental governing equation is the Darcy law. High quality velocity fields are of paramount importance for the successful subsequent coupling with other physical phenomena such as transport. The virtual element method, as solution scheme, is known to be applicable in problems whose discretizations requires cells of general shape, and the mixed formulation is here preferred to obtain accurate velocity fields. To overcome the issues associated to the complex geometries and, at the same time, retaining the quality of the solutions, we present here the virtual element method to solve the Darcy problem, in mixed form, in presence of curved interfaces in two and three dimensions. The numerical scheme is presented in detail explaining the discrete setting with a focus on the treatment of curved interfaces. Examples, inspired from industrial applications, are presented showing the validity of the proposed approach.

Abstract Modern engineering applications require the solution of linear systems of millions or even billions of equations. The solution of the linear system takes most of the simulation for large scale simulations, and represent the bottleneck in developing scientific and technical software. Usually, preconditioned iterative solvers are preferred because of their low memory requirements and they can have a high level of parallelism. Approximate inverses have been proven to be robust and effective preconditioners in several contexts. In this communication, we present an adaptive Factorized Sparse Approximate Inverse (FSAI) preconditioner with a very high level of parallelism in both set-up and application. Its inherent parallelism makes FSAI an ideal candidate for a GPU-accelerated implementation, even if taking advantage of this hardware is not a trivial task, especially in the set-up stage. An extensive numerical experimentation has been performed on industrial underground applications. It is shown that the proposed approach outperforms more traditional preconditioners in challenging underground simulation, greatly reducing time-to-solution.

Abstract Full-physics subsurface simulation models coupled with surface network can be computationally expensive. In this paper, we propose a physics-based subsurface model proxy that significantly reduces the run-time of the coupled model to enable rapid decision-making for reservoir management. In the coupled model the subsurface reservoir simulator generates well inflow performance relationship (IPR) curves which are used by the surface network model to determine well rates that satisfy surface constraints. In the proposed proxy model, the CPU intensive reservoir simulation is replaced with an IPR database constructed from a data pool of one or multiple simulation runs. The IPR database captures well performance that represents subsurface reservoir dynamics. The proxy model can then be used to predict the production performance of new scenarios – for example new drilling sequence – by intelligently looking up the appropriate IPR curves for oil, gas and water phases for each well and solving it with the surface network. All necessary operational events in the surface network and field management logic (such as facility constraints, well conditional shut-in, and group guide rate balancing) for the full-coupled model can be implemented and honored. In the proposed proxy model, while the reservoir simulation component is eliminated for efficiency. The entirety of the surface network model is retained, which offers certain advantages. It is particularly suitable for investigating the impact of different surface operations, such as maintenance schedule and production routing changes, with the aim of minimizing production capacity off-line due to maintenance. Replacing the computationally intensive subsurface simulation with the appropriate IPR significantly improves the run time of the coupled model while preserving the essential physics of the reservoir. The accuracy depends on the difference between the scenarios that the proxy is trained on and the scenarios being evaluated. Initial testing with a complex reservoir with more than 300 wells showed the accuracy of the proxy model to be more than 95%. The computation speedup could be an order of magnitude, depending largely on complexity of the surface network model. Prior work exists in the literature that uses decline curves to replicate subsurface model performance. The use of the multi-phase IPR database and the intelligent lookup mechanism in the proposed method allows it to be more accurate and flexible in handling complexities such as multi-phase flow and interference in the surface network.