Simulating Subsidence from Oil and Gas Extraction
Abaqus FEA can be used to predict ground subsidence and rock deformation of oil and gas reservoirs to enhance safety and operability.
Oil and gas remain primary power sources for both personal and industrial use worldwide. Extraction of these fuel resources from underground reservoirs involves complex geomechanical processes, and can result in subsidence of the ground over a reservoir. Since this occurrence can have an impact on the environment and affect the operability of extraction equipment, it needs to be accurately predicted and kept within safe limits.
During the production, oil and gas pressure is released or, more generally, altered. The fluid pressure distribution in space and time is primarily determined by the permeability and compaction properties of the rock, the locations of the bore holes and the specific extraction schedule and procedures employed. The changes in the pore pressure of the rock, induced by the hydrocarbon extraction, induce changes in the stress state and, as a consequence, deformation of the reservoir rock. This, in turn, leads to modification of the rock’s porosity and permeability, further affecting fluid flow. Simulation of this coupled system requires taking into account both the fluid flow in the pore space and the deformations in the rock (Figure 1).
Abaqus (from SIMULIA, Dassault Systèmes) finite element analysis (FEA) is well suited for modeling the elastic and inelastic deformations in the rock, but it’s not able to calculate the multiphase complexities of the fluid flow regime. Other software tools are able to address the flow phenomena—such as multiple phases, phase changes and miscibility—but they can’t directly address the accompanying deformation and its subsequent impact on flow. To perform a simulation that accurately represents the behavior of the entire system, the two software tools must be combined: Abaqus to compute the deformation in the rock; and a flow simulation tool (such as ECLIPSE from Schlumberger) to calculate flow and the pore-pressure depletion history.
Studies that used such a coupled approach required considerable manual effort to create and modify the finite element mesh and transfer the pore pressure data from ECLIPSE to Abaqus/Standard. A more recent study demonstrated how capabilities that are now available in Abaqus can be used to automate several of these tasks.
For the current analysis, the PUNQ reservoir was studied using an automatic workflow approach to carry out field-scale geomechanical modelling.
As a first step, the ECLIPSE flow simulation grid and result data were used to generate an Abaqus/CAE output database using a translator, obtaining an equivalent model in terms of grid and properties needed for the geomechanical analysis. Initial values of porosity as well as the pore pressure depletion history were automatically assigned to the elements generated in the output database.
Although multiple fluids exist within each flow simulation grid cell, Abaqus is able to represent only a single fluid within each finite element. To select the appropriate fluid for each element, the translator looked at the ECLIPSE data and automatically designated regions that were principally gas-, oil- or water-bearing for each reservoir layer, while also deriving specific weight values for the various fluids. These results were written to the same database, which was then imported into Abaqus/CAE for further analysis (Figure 2).
Refining the model, a script in Abaqus was used to manage a sequence of operations: first, the mesh was modified to merge extremely thin layers into a single layer; then density, specific weight, material properties, initial and boundary conditions were specified. Since the flow simulation grid mesh was different than the mesh to be used for the geomechanics analysis (all cell-centered properties), variables (including pore pressure and initial void ratios) needed to be accurately transferred from one mesh to the other. Mapping functionalities in Abaqus were utilized to transfer the data.
Layers of elements were next added to the model to represent the under-, side- and over-burden regions, the grid for these regions not being present in the original flow simulation model (Figure 3). The burden regions were added such that the top surface of the model was at ground level, the bottom was at a depth of 5 km, and the final surface dimensions were 15.0 by13.5 km.
Final model set-up included the following steps: Initial stresses were represented as a piece-wise distribution through the model depth; initial pore pressures and void ratios were read directly from the output database created by the translator; displacements normal to the side and bottom boundary surfaces were specified as zero; material properties were assigned directly in CAE using an “ad hoc” developed procedure; pore pressure boundary conditions (determined through a sub-modeling technique) were driven by values computed from the flow simulation and available in the output database and, since pore pressures were specified at all nodes having pore pressure degrees of freedom, the permeability values were deemed inconsequential.
The coupled simulation involved two distinct analyses. The first—an elastic geostatic analysis—was performed to obtain the vertical stress distribution. For this calculation, the porous reservoir rock was modeled as porous elastic, while the non-porous regions (the under-, side- and over-burden regions) were modeled as linear elastic. The software iterated until the displacements obtained for the applied gravity load were nearly zero. The result was a stress distribution that equilibrated the applied loading and boundary conditions.
For the second analysis—in which the reservoir rock was modeled as elastoplastic—vertical stress values from the first analysis were used to compute the stress-dependent compressibility values of the reservoir rock, which were then utilized to update the linear elastic properties of the non-porous regions. Next, a steady-state coupled pore fluid flow displacement calculation was performed for the time increments representing the pore-pressure depletion data sets obtained from the flow simulation. This calculation yielded displacement and plastic strain values for the reservoir and surrounding burden regions.
A selection of the results achievable with the analyses is presented in the following, chosen as the most interesting with relation to the subsidence study:
- The downwards settlement of top surface of the reservoir (Figure 4a);
- The plastic strain magnitude: its increase with time indicates that the porous region is compacting as fluids is extracted (Figure 4b);
- The maximum subsidence value: it occurred immediately above the reservoir and decreased rapidly as distance from this location increased (Figure 4c).
The results of the analysis showed a strong correlation to those of the first analysis, which used manual methods to transfer data and modify the mesh.
The methodologies described here can be used to predict reservoir compaction and surface subsidence of the ground. The geomechanics model can serve as a global model in the creation of sub-models for smaller scale applications, such as bore-hole stability. Additionally, the study validated that its methods can be used as one of the components of an iteratively coupled geomechanics simulation.