A consistent framework for coupling modern reservoir flow simulator with mechanics

This work presents a numerical framework, for the analysis of coupled fluid and geomechanical processes related to subsurface operations, built upon ABAQUS of SIMULIA, an industrial-standard software for mechanical analysis, and Echelon, a GPU-based reservoir simulator. By using porosity updates, the proposed framework can address both the impact of pore pressure changes on rock compaction as well as the the feedback of skeleton deformations on multiphase flow. Our solution is validated by means of various test cases. The results agree with analytical and reference solutions, showing that the proposed implementation of the iterative algorithm correctly accounts for strong coupling between mechanics and flow in porous media.


Introduction
The depletion of hydrocarbon reservoirs induces variations in time and space of reservoir pressure and saturation. In turn, changes in reservoir hydraulic properties may cause, following Terzaghi's stress principle, a modification of the stress state in and around the reservoir. Stress changes can then influence the reservoir fluid parameters and alter the production scenario with the possibility of well mechanical failure, induced fracturing, fault reactivation and ground surface deformations. Nonetheless, the standard practice in the Oil&Gas industry is to assume that porosity -and permeability -are function of fluid pressure only and not of the stress state. This may not capture the coupling between flow and geomechanics, specially for weak hydrocarbon reservoirs, where the complex stress state and loading conditions give rise to complicated constitutive behaviour and stress-dependent reservoir properties [10,11,8].
Coupled hydro-mechanical analysis can be performed using either fully coupled or sequential methods. Unknowns for flow and mechanics are either solved simultaneously within a timestep (fully coupled or monolithic), or sequentially by splitting the problem into sub-problems. Fully coupled methods are attractive for general purpose simulations due to its unconditional stability [14]. However, they are computationally expensive and require the development of a unified flow-geomechanics simulation framework inside a single code.
On the other hand, sequential methods represent a viable compromise between accuracy and computational efficiency. They have the advantage to allow the use of separate packages for each physical sub-problem that can be solved with dedicated numerical methods and optimum spatial discretization in the different spatial domains. For instance, conventional reservoir simulators are based on first order, cell centered, Finite Volume (FV) schemes, while stress equilibrium is usually computed by means of Finite Element (FE). The simplest sequential scheme is the IOP Conf. Series: Earth and Environmental Science 833 (2021) 012113 IOP Publishing doi:10.1088/1755-1315/833/1/012113 2 one-way coupling: the reservoir simulator is run for a prescribed time period and afterwards the calculated pressure is provided as input for the finite element analysis to evaluate induced changes in the stress state, rock compaction and subsidence. Notably, one-way coupling is easy to implement and it is effective for a number of problems, particularly subsidence studies for gas reservoirs. However, in case of strongly coupled systems as in stress-dependent reservoirs, oneway coupling is known to give erroneous results [11,13,9]. A more accurate approach is given by the two-way coupling method where information between the flow and geomechanic simulators is transferred back and forth until convergence of both problems is reached. Convergence rate and stability of the iterative method strongly depend on the sequential scheme: the fixed-stress split has been proved to be unconditionally stable and to converge to the fully coupled solution [4,1].
In this work, we propose the implementation of a consistent iterative framework -based on the fixed-stress split proposed in [4] -to couple industrial-grade reservoir simulators with geomechanical solvers. The framework couples ABAQUS, distributed by Simulia, and Echelon, a GPU-based reservoir simulator [6] jointly developed by Eni and Stoneridge Technology [2]. We developed a communication interface that allows for accounting for both the influence of pore pressure on the stress equilibrium and that of deformations on multiphase flows. Throughout socket technology, we use the Co-Simulation Engine (CSE) developed by Simulia for the communication with ABAQUS and standard Trasmission Control Protocol and Internet Protocol (TCP/IP) for the communication with Echelon.
In the next section we briefly summarize the mathematical models of the different physical processes involved (flow and mechanics), then, in Sec. 3, we describe the implementation of the fixed-stress scheme and the solution strategy and finally, in Sec. 4, we introduce a set of numerical examples to validate the proposed coupling framework.

Mathematical model
In this section, we summarize the mathematical models that describe the mass balance for multiple fluid phases in porous media, the momentum balance for the rock matrix and the corresponding coupling parameters.

Mass conservation equations for fluid phases
The black-oil model is the most used formulation in the Oil&Gas community [5]. It gives the possibility to model mass transfer between oil and gas phases at reservoir conditions assuming conservation of stock tank components with an immiscible water phase. Mass conservation equation reads as where α labels water, oil and gas phases, φ is the porosity, K is the absolute permeability tensor, S α , ρ α , k r,α and µ α are the phase saturation, density, relative permeability and viscosity, respectively, and Φ α = p α − ρ α gz is the phase potential, which corrects the phase pressure p α with a gravity contribution. Fluids fully saturate the pore space, S w + S o + S g = 1, and capillary pressure function relates phase pressure to saturation, p o −p w = p c,ow (S w ) and p g −p o = p c,og (S g ). Flow equations are solved using a first order FV discretization, where primary variables are cell pressure and water/gas saturations.

Momentum conservation equations for solid
Assuming quasi-static conditions, the conservation of linear momentum for the fully saturated rock matrix is: where σ is total Cauchy stress tensor and ρ b = ρ s (1 − φ) + φρ f is the overall (bulk) mass density. Here, ρ s is the grain density and ρ f = α S α ρ α is the average fluid density. Following Terzaghi's stress principle and Biot's theory [12], the total stress can be decomposed into an effective stress and a fluid pressure contribution: Here, C is the fourth-order elasticity tensor, = ∇u + ∇ T u /2 the second-order strain tensor, u the displacement vector, b the Biot coefficient and, P the saturation average fluid pressure. The convention on the stress and pressure sign is negative for compression. For isostropic, linear elastic materials, the stress/strain relationship reads as Assuming that fluid pressure is known, Eqs. 2 and 3 are discretized by FE, where the unknowns are the nodal displacements.

Coupling parameters
Multiphase flow and geomechanics are coupled through the following relationship between porosity and stress/strain [12,4]: where v = tr( ) is the volumtric strain, K s is the grain bulk modulus and φ 0 , P 0 and v,0 denote reference porosity, pore pressure and volumetric strain, respectively. According to Eq. 4, strain (or effective stress) variations cause changes in the pore volume, which accordingly influence the pressure of pore filling fluids, while the ratio of fluid volume change induced by bulk volume changes is expressed by the Biot's coefficient.

Implementation of Sequential Coupling
Sequential coupling envisages iterative solution of the two sub-problems, with cell pressure (p) and saturations (S w/g ) as unknowns of the flow problem, while nodal displacement (u) as unknowns for mechanics. Solution is then achieved as follow: for every time-step n and coupling iteration k, the flow problem is solved first using the volumetric strain from the previous iteration to compute the porosity correction. In order to decouple the flow problem from mechanics, we employ the fixed-stress split, which provides an unconditionally stable sequential formulation [4,1]. The volumetric total stress, σ v , is constant between iteration k and k + 1, i.e.
where we have assumed negligible capillary pressure so that P = p o . Eq. 4 at iteration k + 1 is Inserting Eq. 5 into Eq. 6, we obtain the final form of the porosity correction: The sequential algorithm is implemented by means of a communication interface (C++) which connects the two simulators, namely ABAQUS and Echelon, and defines when convergence is achieved. The implementation is based on standard TCP/IP technology, where sockets are used to minimize simulator code modifications. After an initial equilibration step where constant rock properties are assumed, the controller drives the simulation from time t n to time t n+1 as follows: where ε p,u are user-defined tolerances, the simulation is advanced to the next time-step, otherwise a new iteration is performed.
Since ABAQUS computes the displacements at the nodes of each finite element, while pressure is assigned at the center of each cell, an interpolation step is necessary. For each node of the FE grid, a cell-to-node map is constructed during pre-processing and then applied during solution to interpolate pressure from cell center to cell nodes.

Numerical examples
In this section, we validate the proposed implementation of the coupling algorithm and demonstrate its applicability to multi-physics problems of practical interest. We use three test cases for this purpose. In the first two cases we compare our numerical results with available analytical or numerical benchmarks for single-phase problems. The third case demonstrates that the sequential approach can accurately solve a realistic multi-phase scenario of injection and production. In this case, the validation is based on a fully coupled solution provided by AD-GPRS, an academic code developed by Standford University [3]. For the purpose of validation, Echelon and ABAQUS share the same reservoir grid in all the proposed test cases, which makes the construction of the cell-to-node map for pressure interpolation almost straightforward. Note however, that the current implementation is not limited to the same grid for flow and mechanics as long as a suitable interpolation map is provided. Moreover, the geomechanical grid can be extended to include overburden, underburden and sidesburden by assuming that the additional finite elements are either non-porous or at constant (hydrostatic) pore pressure.

Mandel's problem
Mandel's problem [15] consists of an homogeneous poroelastic medium sandwiched between two rigid, frictionless, and impermeable plates, both loaded by a constant perpendicular force. The slab is infinitely long in the y-direction with a 2a × 2b wide cross-section. The lateral surfaces are traction-free and drained. The overburden load is suddenly applied at t = 0, causing both an instantaneous overpressure and settlement. Due to the symmetry of the problem, only one quarter of the domain is considered. The medium is discretized with 30 × 10 gridblocks with ∆x = ∆z = 1 m. We use a uniform initial pressure p 0 = 2.125 M P a and we characterized mechanics and flow properties as follows:

Dean's problem
This validation problem was initially analyzed by [10] and afterwards presented by [13]. It consists of a single-phase, soft reservoir that is contained in a stiff non-pay region. The Young moduli of the reservoir and non-pay region differ by two orders of magnitude, which enhances the coupling effects and gives rise to geomechanical responses at the boundary of the reservoir that can not be capture with uncoupled simulations (Mandel-Cryer effect). The reservoir grid consists of 11 × 11 × 5 gridblocks included at the center of a larger mechanical grid (21 × 21 × 12 gridblocks), with five layers of overburden, two layers of underburden and five layers on each side of the reservoir of sideburden (see Fig. 2a). A single vertical well is completed in the center of the reservoir in all five layers and produces at constant rate of 7950 sm 3 /d. Water production is simulated for 3000 days, with a time-step size of 20 days for the first 400 days, followed by timesteps of 200 days till the end of the simulation. A detailed description of the problem, including grid size, fluid and rock properties can be found in [13]. Figs. 2b and 2c show a comparison of the simulated pressure at the boundary of the reservoir in cell (6,11,6) and the subsidence at the surface together with the reference solutions. For comparison purposes, the one-way-approach results, obtained by limiting the number of coupling iteration to one, are also shown in Fig. 2. The one-way approach underestimates the pressure drop in the reservoir and the subsidence at the surface. As observed in the Mandel's problem, the proposed coupled method captures the non-monotonic pressure response due to fast stress equilibration at the boundary of the reservoir confirming the accuracy of the implementation.
We investigate the convergence of the coupling algorithm at different coupling strengths. As outlined in [4], the coupling strengths is defined as the ratio of the bulk stiffness of the fluid and solid skeleton τ = b 2 M/K d , where M −1 = φc f + (b − φ)c s is the inverse of the storativity coefficient and accounts for both the fluid compressibility, c f , and grain compressibility, c s = K −1 s . Large values of τ indicate strongly coupled systems characterized by soft material and low compressibilities, for which strong geomechanical influence on the fluid flow are expected. The inset of Fig. 2d shows the pressure behaviour at the boundary of the reservoir for different values of τ . At increasing coupling strengths, i.e. decreasing grain compressibility in this case, pressure depletion due to production is faster and the difference with respect to the uncoupled case increases. Note that the uncoupled case, as in traditional  Figure 2: Dean's problem: (a) geometry, (b) pore pressure in cell (6,11,6) and (c) subsidence as a function of time for coupled (two-way) and uncoupled (one-way) simulations. Numerical results are shown with circles, while dashed lines represent the reference solutions [13]. (d) Number of coupling iterations per time step and (inset) reservoir pressure in cell (6,11,6) at different coupling strengths.
reservoir simulations, is independent of the grain compressibility. At the same time, as shown in Fig. 2d, at increasing coupling strengths more iterations are required to solve the coupled system and thus more computational efforts.

Synthetic two-phase problem
In the last numerical example we consider a synthetic two-phase problem. The reservoir consists of an under-saturated oil accumulation in a four-way anticline structure surrounded by a stiff non-pay region, see Fig. 3. A corner point geometry grid, see [7], consisting of 28 × 27 Water compressibility is 4 10 −4 MPa −1 , ρ w = 1000 kg/m 3 and µ w = 10 −3 kg/ms. Before production, reservoir is in equilibrium with a fluid pressure of 30 kg MPa at oil-water contact. We use quadratic relative permeabilities with a critical water saturation and critical oil to water saturation of 0.2.
The reservoir is exploited by 3 wells, 2 producers and 1 injector, with a primary depletion (4500 days) followed by 4500 days of water injection. Wells produce with a maximum liquid rate limit of 7950 sm 3 /d and a minimum bottom-hole pressure of 20 MPa, while the injector operates at a maximum bottom-hole pressure of 350 MPa. Figure 4: Results of the synthetic two-phase problem for the coupled (two-way) and uncoupled (one-way) cases: (left) average hydrocarbon pressure and (right) pressure in cell (1, 14, 1). Oneway and fully-coupled solutions obtained with AD-GPRS are superimposed as references.
Average pressure and pressure in cell (1, 14, 1) along of time are shown in Fig. 4 together with AD-GPRS fully coupled results. As in the Dean's problem, one-way coupling approach underestimates the pressure drop in the reservoir and can not capture Mandel-Cryer effect at the reservoir boundary. Moreover, the comparison with AD-GPRS shows that the results obtained with our sequential coupling approach converge to the fully-coupled solution. The case was run using a single NVIDIA Tesla-V100 GPU for Echelon and 8 processors (Intel Xeon Silver 4116, 10GHz ) for ABAQUS. In this case, the computational efforts required for solving the fluid dynamic problem and for communication are negligible with respect to those required for the solution of the mechanical problem (more than 99% of the total simulation time). Further code optimizations are left for future studies.

Conclusion
We presented a consistent framework, based on the fixed-stress sequential approach, for coupling a state-of-the-art reservoir simulator with a geomechanical solver. The coupling framework is built upon ABAQUS and Echelon, a GPU-based reservoir simulator. The implementation, which accounts for both the impact of pore pressure changes on rock stability and deformation as well as the feedback of rock compaction on the fluid behaviour, hinges on a communication interface developed to exchange pressure and stress at each coupling iteration between the two simulators. The communication package is modular, which provides a viable solution to integrate a generic geomechanical software with the reservoir simulator.
The coupling algorithm is validated against several reference benchmarks. More specifically, it is shown that our solution is suitable to correctly model the Mandel-Cryer effect, a nonmonotonic pressure response to fast stress re-equilibration, in problems where analytical or reference solutions are available and in realistic field-like problems. Finally, we have shown that our implementation converges to the fully coupled solutions.