In this research we consider developing a reservoir simulator capable of simulating complex coupled poromechanical processes on massively parallel computers. A variety of problems arising from petroleum and environmental engineering inherently necessitate the understanding of interactions between fluid flow and solid mechanics. Examples in petroleum engineering include reservoir compaction, wellbore collapse, sand production, and hydraulic fracturing. In environmental engineering, surface subsidence, carbon sequestration, and waste disposal are also coupled poromechanical processes. These economically and environmentally important problems motivate the active pursuit of robust, efficient, and accurate simulation tools for coupled poromechanical problems. Three coupling approaches are currently employed in the reservoir simulation community to solve the poromechanics system, namely, the fully implicit coupling (FIM), the explicit coupling, and the iterative coupling. The choice of the coupling scheme significantly affects the efficiency of the simulator and the accuracy of the solution. We adopt the fixed-stress iterative coupling scheme to solve the coupled system due to its advantages over the other two. Unlike the explicit coupling, the fixed-stress split has been theoretically proven to converge to the FIM for linear poroelasticity model. In addition, it is more efficient and easier to implement than the FIM. Our computational results indicate that this approach is also valid for multiphase flow. We discretize the quasi-static linear elasticity model for geomechanics in space using the continuous Galerkin (CG) finite element method (FEM) on general hexahedral grids. Fluid flow models are discretized by locally mass conservative schemes, specifically, the mixed finite element method (MFE) for the equation of state compositional flow on Cartesian grids and the multipoint flux mixed finite element method (MFMFE) for the single phase and two-phase flows on general hexahedral grids. While both the MFE and the MFMFE generate cell-centered stencils for pressure, the MFMFE has advantages in handling full tensor permeabilities and general geometry and boundary conditions. The MFMFE also obtains accurate fluxes at cell interfaces. These characteristics enable the simulation of more practical problems. For many reservoir simulation applications, for instance, the carbon sequestration simulation, we need to account for thermal effects on the compositional flow phase behavior and the solid structure stress evolution. We explicitly couple the poromechanics equations to a simplified energy conservation equation. A time-split scheme is used to solve heat convection and conduction successively. For the convection equation, a higher order Godunov method is employed to capture the sharp temperature front; for the conduction equation, the MFE is utilized. Simulations of coupled poromechanical or thermoporomechanical processes in field scales with high resolution usually require parallel computing capabilities. The flow models, the geomechanics model, and the thermodynamics model are modularized in the Integrated Parallel Accurate Reservoir Simulator (IPARS) which has been developed at the Center for Subsurface Modeling at the University of Texas at Austin. The IPARS framework handles structured (logically rectangular) grids and was originally designed for element-based data communication, such as the pressure data in the flow models. To parallelize the node-based geomechanics model, we enhance the capabilities of the IPARS framework for node-based data communication. Because the geomechanics linear system is more costly to solve than those of flow and thermodynamics models, the performance of linear solvers for the geomechanics model largely dictates the speed and scalability of the coupled simulator. We use the generalized minimal residual (GMRES) solver with the BoomerAMG preconditioner from the hypre library and the geometric multigrid (GMG) solver from the UG4 software toolbox to solve the geomechanics linear system. Additionally, the multilevel k-way mesh partitioning algorithm from METIS is used to generate high quality mesh partitioning to improve solver performance. Numerical examples of coupled poromechanics and thermoporomechanics simulations are presented to show the capabilities of the coupled simulator in solving practical problems accurately and efficiently. These examples include a real carbon sequestration field case with stress-dependent permeability, a synthetic thermoporoelastic reservoir simulation, poroelasticity simulations on highly distorted hexahedral grids, and parallel scalability tests on a massively parallel computer. / text

There are two big challenges which restrict the extensive application of fully coupled geomechanics-reservoir modeling. The first challenge is computational effort. Consider a 3-D simulation combining pressure and heat diffusion, elastoplastic mechanical response, and saturation changes; each node has at least 5 degrees of freedom, each leading to a separate equation. Furthermore, regions of large p, T and σ′ gradients require small-scale discretization for accurate solutions, greatly increasing the number of equations. When the rock mass surrounding the reservoir region is included, it is represented by many elements or nodes. These factors mean that accurate analysis of realistic 3-D problems is challenging, and will so remain as we seek to solve larger and larger coupled problems involving nonlinear responses. To overcome the first challenge, the displacement discontinuity method is introduced wherein a large-scale 3-D case is divided into a reservoir region where Δp, ΔT and non-linear effects are critical and analyzed using FEM, and an outside region in which the reservoir is encased where Δp and ΔT effects are inconsequential and the rock may be treated as elastic, analyzed with a 3D displacement discontinuity formulation. This scheme leads to a tremendous reduction in the degrees of freedom, yet allows for reasonably rigorous incorporation of the reactions of the surrounding rock. The second challenge arises from some forms of numerical instability. There are actually two types of sharp gradients implied in the transient advection-diffusion problem: one is caused by the high Peclet numbers, the other by the sharp gradient which appears during the small time steps due to the transient solution. The way to eliminate the spurious oscillations is different when the sharp gradients are induced by the transient evolution than when they are produced by the advective terms, and existing literature focuses mainly on eliminating the spurious spatial temperature oscillations caused by advection-dominated flow. To overcome the second challenge, numerical instability sources are addressed by introducing a new stabilized finite element method, the subgrid scale/gradient subgrid scale (SGS/GSGS) method.

There are two big challenges which restrict the extensive application of fully coupled geomechanics-reservoir modeling. The first challenge is computational effort. Consider a 3-D simulation combining pressure and heat diffusion, elastoplastic mechanical response, and saturation changes; each node has at least 5 degrees of freedom, each leading to a separate equation. Furthermore, regions of large p, T and σ′ gradients require small-scale discretization for accurate solutions, greatly increasing the number of equations. When the rock mass surrounding the reservoir region is included, it is represented by many elements or nodes. These factors mean that accurate analysis of realistic 3-D problems is challenging, and will so remain as we seek to solve larger and larger coupled problems involving nonlinear responses. To overcome the first challenge, the displacement discontinuity method is introduced wherein a large-scale 3-D case is divided into a reservoir region where Δp, ΔT and non-linear effects are critical and analyzed using FEM, and an outside region in which the reservoir is encased where Δp and ΔT effects are inconsequential and the rock may be treated as elastic, analyzed with a 3D displacement discontinuity formulation. This scheme leads to a tremendous reduction in the degrees of freedom, yet allows for reasonably rigorous incorporation of the reactions of the surrounding rock. The second challenge arises from some forms of numerical instability. There are actually two types of sharp gradients implied in the transient advection-diffusion problem: one is caused by the high Peclet numbers, the other by the sharp gradient which appears during the small time steps due to the transient solution. The way to eliminate the spurious oscillations is different when the sharp gradients are induced by the transient evolution than when they are produced by the advective terms, and existing literature focuses mainly on eliminating the spurious spatial temperature oscillations caused by advection-dominated flow. To overcome the second challenge, numerical instability sources are addressed by introducing a new stabilized finite element method, the subgrid scale/gradient subgrid scale (SGS/GSGS) method.