SPE-182721-MS A Comprehensive Computational Framework for Long Term Forecasting on Serious Wellbore Damage and Ground Surface Subsidence Induced by Oil and Gas Reservoir Depletion Ruijie Liu and Zhijun Liu, Department of Mechanical Engineering, University of Texas at San Antonio Copyright 2017, Society of Petroleum Engineers This paper was prepared for presentation at the SPE Reservoir Simulation Conference held in Montgomery, TX, USA, 20–22 February 2017. This paper was selected for presentation by an SPE program committee following review of information contained in an abstract submitted by the author(s). Contents of the paper have not been reviewed by the Society of Petroleum Engineers and are subject to correction by the author(s). The material does not necessarily reflect any position of the Society of Petroleum Engineers, its officers, or members. Electronic reproduction, distribution, or storage of any part of this paper without the written consent of the Society of Petroleum Engineers is prohibited. Permission to reproduce in print is restricted to an abstract of not more than 300 words; illustrations may not be copied. The abstract must contain conspicuous acknowledgment of SPE copyright. Abstract Depletion of oil and gas reservoirs results in ground surface subsidence and wellbore damage. Accurate prediction of subsidence and wellbore damage can provide a strong support to reservoir management and community precaution on subsidence impact. For reservoirs with soft rock formations, depletioninduced reservoir compaction results in serious damage and permanent deformation of rock formation. Cap plasticity models are crucial to achieve a more reliable prediction on the compaction of reservoirs with soft rocks. In this study, we develop a fully coupled and fully implicit finite element framework along with cap plasticity models for modeling the compaction of reservoirs with soft formations undergoing finite deformation. Creep models are also formulated with cap plasticity models for predicting the longterm effect of reservoir compaction. A material integrator for stress return mapping for both plasticity and creep models is consistently formulated to achieve fast convergent rates. A reservoir compaction model with three horizontal production wells and a wellbore damage model are solved using the proposed computational framework as well as cap and creep models with parameters obtained from field tests. We have demonstrated: 1) a good performance of the proposed computational framework in modeling the compaction of soft rock reservoirs with complex payzone geometries; and 2) the predicted subsidence with cap plasticity models could be seven times greater than linear elastic or plasticity models without caps. Introduction Depletion of oil and gas reservoirs results in rock compaction, wellbore instability, and ground surface subsidence. These phenomena have been making persistent unfortunate impacts on oil industries and societies. Assessment and prediction on these reservoir depletion induced events are crucial for improving reservoir management and preventing potential man-made disasters. For reservoirs with soft rocks undergoing serious compaction, formation damage and failure progressively develop in pay zones. Numerical simulation and analysis on subsidence (Rattia et al. 1981, Bruno 1992) and wellbore damage (Vudovich et al. 1989) have been extensively studied for Ekofisk (Baade et al. 1988, Chin et al. 1993), North sea (Keszthelyi et al. 2016) and the Belridge diatomite field (De Rouffignac et al. 1995, Hansen et al. 1995). 2 SPE-182721-MS Popular computational frameworks (Settari et al. 1998, Settari et al. 2001, Gai et al. 2003) in modeling reservoir subsidence are based on linear elasticity or Drucker-Prager plasticity models within the framework of small deformation. Different coupling methods for reservoir modeling are studied in (Dean et al. 2006), including fully coupled method, explicit coupling method, and iterative coupling methods. For reservoirs with soft rock formations, serious subsidence and wellbore failure are mainly caused by large and permanent deformation in terms of formation compaction. Geomechanics cap models have been developed (Fossum et al. 2000, Settari et al. 2001, Foster et al. 2005, Liu et al. 2009, Sun et al. 2014, Motamedi et al. 2015) for modeling this compaction dominant deformation for geomaterials. In (Settari et al. 2001), a Drucker-Prager plasticity model with cap is implemented through designing a transitional zone that smoothly connects the shear envelope and the compaction zone in the yield function. Another important cap model was proposed by Pelessone (Pelessone 1989) by introducing a Heaviside function to bend the shear envelope toward the elliptic yield surface to construct a smooth yield surface for cap models. The Pelessone smooth cap model has been further extended by the incorporation of the stress triaxiality and the Lode angle and also strong hardening functions for the compaction cap portion (Foster et al. 2005). Consistent finite element formulations with small deformation have been constructed in (Fossum et al. 2000, Foster et al. 2005) for the Pelessone-based smooth cap models. To achieve quadratic convergent rates, completely implicit material integrators for stress return mapping have been rigorously constructed in these efforts by performing consistent linearization. Good performance was demonstrated within pure solid field using simple problems with a single element or problems with shear dominant behaviors. More complex, compaction-dominant, and field-like numerical examples with sharp cap models successfully solved by these implicit formations have not been reported yet. A coupled fluid flow and geomechanics with the Pelessone smooth cap model was formulated in (Sun et al. 2014) to model hydro-mechanical responses of strip and circular punching loading for geomaterials. In (Sun et al. 2014), to achieve convergent solutions an semi-implicit formulation was adopted for material integrators for cap models. However, this semi-implicit formation is not the preferred computational framework because stresses are not in exact stress spaces and results may have large errors for solutions at long time. Few numerical examples from current literatures presented large and serious subsidence predictions that were dominated by rock formation compaction failure. On the other hand, for serious compaction cases large deformation occurs and computational frameworks for small deformation are not appropriate. In this work, to accurately model the compaction of reservoirs with soft rocks, a fully coupled and fully implicit finite element framework along with the Pelessone smooth cap model combined with a creep model has been developed for solving problems with finite deformation. The outline of the paper is as follow. In the next section, governing equations for the coupled fluid flow and geomechanics problems as well as the corresponding boundary and initial conditions are presented. In the third section, we briefly review the Pelessone cap plasticity theory and also introduce a creep model. Finite element formulation and implementation for cap and creep models are summarized in the fourth section. In the fifth section, two numerical examples, a reservoir model with three horizontal production wells and a wellbore model are solved using our proposed computational framework. Conclusions are summarized in the last section. Governing Equations Governing equations for coupled fluid flow and geomechanics problems can be summarized as follows: (1) with boundary conditions: SPE-182721-MS 3 (2) and initial conditions: (3) in Eqs. (1)–(3), Ω is the domain, u the displacement, P the pore pressure, f the body force, σ the total stress, σ' the effective stress for solid skeletons. The total stress of porous media filled with fluids is decomposed into pore pressure for fluids and effective stress for solid skeletons as follows: (4) where I is the second-order identity tensor. α is the Biot's constant scaler parameter. In the flow equation in Eq. (1), kp is the permeability of the coupling system, μ the viscosity of fluids, and M the Biot's modulus (Biot 1941). M can be further computed as follows: (5) where Ks is the bulk modulus of solid grains, Kf the bulk modulus of fluids, and ϕf the porosity. Prescribed displacement boundary condition and traction boundary condition for the solid field are presented in the first and second equations in Eq. (2). The last two equations in Eq. (2) summarize the prescribed pressure and flux boundary conditions for the fluid flow field. In Eq. (3), initial total stresses for porous media and initial pore pressure for fluid flow field are prescribed. Cap Plasticity and Creep Models Traditional Drucker-Prager plasticity model is suitable for capturing the behavior of geomaterials dominated by shear plastic deformation. For geomaterials subject to compaction, cap plasticity models are crucial. As non-smoothed cap models (Sandler et al. 1976) have corners in yield surfaces which bring singularities at the intersection of shear and cap part of yield functions, and may often result in non-convergent solutions, smooth cap models (Desai et al. 1986, Pelessone 1989, Schwer et al. 1994) have been developed to avoid these issues. In this study, we adopt the Pelessone smooth cap model (Pelessone 1989, Foster et al. 2005). In addition, creeping effect is important for modeling reservoir compaction for long time. We also introduce a creep model in our proposed computational framework. Yield and Flow Functions The yield function (Foster et al. 2005) is (6) where 4 SPE-182721-MS (7) (8) (9) (10) In Eqs. (6)–(10), I1, J2 and J3 are the first invariant of stress tensor, the second, and the third invariants of deviatoric stress tensor, respectively. X is the intersection of yield function with I1 axis in the compression portion in meridional space, k the transition between shearing and compaction portions of the yield function. Af, Bf and αf are shearing yield related parameters, σ0 the cohesion related shearing yield parameter. ψ is the is the Heaviside ratio of the strength in tension and compression. Rf determines the shape of the cap. function. An example of the schematic of the yield function is shown in Fig. 1. Figure 1—Schematic of the yield function. (a) 3-d plot of yield function in principle stress space; (b) plot of yield function in meridional coordinate system. The expression of flow function is basically the same with the yield function, except the parameters with subscript "f" in Eqs. (6)–(10) are replaced with their counterparts with subscript "g". Hardening The hardening function of σ0 is described by (11) where ∊eq, in is the equivalent inelastic strain. The hardening of X is expressed by the following equation: (12) SPE-182721-MS 5 where is the inelastic volumetric strain. The relation between cap hardening parameter X and inelastic volumetric strain is: (13) Creeping A creep model is incorporated in this study. We adopt the creep model in (Rinehart et al. 2014) as follows: (14) are two scalars and they are equivalent creeping strain and stress. and are, in general, where and measured from simple experimental tests and their relations with the general three dimensional stress and creep strain by the creep power are as follows: (15) where ∊c is the creep strain tensor. The creep rate is obtained as follows: (16) where the same plastic potential function is assumed for creep potential. Finite Element Formulation Finite element formulations for coupled fluid flow and geomechanics problems governed by equations in previous sections can be found in (Lewis et al. 1987, Liu et al. 2009, Liu et al. 2009). The algebraic equations are summarized below: (17) where Kuu is the stiffness matrix corresponding to the solid skeleton, which is the sum of material stiffness, geometric stiffness and loading stiffness in finite deformation problems; Kup corresponds to the coupling term. Kpp is the stiffness of the fluid flow. Ru is the residue corresponding to the equilibrium equation, which is caused by both stress and external loading, Rp the residue corresponding to the fluid flow. Subsidence in many reservoirs with soft rocks could be up to more than several ten ft. and volume strain may be up to 20%. In this study, we have formulated the coupled problem within finite deformation frameworks. The co-rotational framework (Oden 1972) is adopted to implement hypoelastoplasticity model. Eq. (16) is highly nonlinear due to plasticity models. To achieve fast convergent rates, the formulation is required to consistently linearize Eq. (16). A key linearization is to construct material integrator to perform stress return mapping. The material integrators provide: 1) updated stress that is on the yield surface if materials yield; and 2) tangent stiffness for the construction of global stiffness matrix. Consistently formulated material integrators for the Pelessone cap plasticity model can be found in (Foster et al. 2005). We have reformulated these formulations for cap plasticity models to take into account creeping behaviors of rock formations. The key procedure for our new formulation for both cap plasticity and creep models is to solve a pure creeping problem first. Material yield condition will be checked after solutions converge for creeping problem. If materials yield, the problem will be re-solved again for combined plasticity and creep 6 SPE-182721-MS models. A computer code "Integrated Geomechanics Analysis of Reservoir (IGAR)" has been developed in the University of Texas at San Antonio that focuses on simulating reservoir compactions. We have performed a consistent formulation in this code for coupled fluid flow and geomechanics problems with both cap and creep models by consistently linearizing Eq. (16). Quadratic convergent rates for caps with sharp transition zones are achieved in this code. Numerical Examples In this section, we present two numerical examples to demonstrate the procedure and the performance of our computational framework for modeling reservoir subsidence and wellbore damage with experiment-based data. First, a reservoir model with complex payzone geometry, sharp cap models, and multiple horizontal wells is solved using our IGAR code. We also solve a wellbore deformation problem using IGAR code in which rock formation is in compaction. Reservoir Subsidence Fig. 2 presents a plain strain reservoir model having three geology layers with 5000 ft. in horizontal direction and 600 ft. in vertical direction. In the model, the top layer represents an impermeable caprock. The payzone in the middle layer has 200 ft. in thickness and is composed of diatomite which are typical formations for the Belridge diatomite field. The bottom layer is an impermeable underburden layer. The sideburdens are constrained in horizontal directions for normal displacements. The unerburden is constrained on the vertical direction for displacements. Pressure loading with 1000 psi is applied on the top overburden of the model. Three horizontal wells are distributed into the payzone. The wellbore pressure around these three wells are assumed as p = 400 psi. For the payzone layer, we set up initial pore pressure as p0 = 900 psi and initial stresses as σx0 = σz0=-5 psi, σy0=-80 psi. Linear elasticity model is assumed for the overburden and underburden. The material parameters including Young's modulus, Poisson ratio, and the Drucker-Prager plasticity with cap model data are listed in Table 1 where data is based on experimental tests and collected from (Fredrich et al. 1996, Rinehart et al. 2014). In the model, a total of 952 linear hexahedrons for both solid and fluid fields are meshed, in which relative fine meshes are used near the wellbores and the center of the reservoir. Figure 2—Reservoir model SPE-182721-MS 7 Table 1—MATERIAL PARAMETERS FOR RESERVOIR MODEL Parameter Description Value Unit E1 Young's modulus of top layer 6.8 × 106 psi E2 Young's modulus of pay zone 6.0 × 104 psi E3 Young's modulus of bottom layer 3.6 × 107 psi ν1 Poisson's ratio of top layer 0.22 Dimensionless ν2 Poisson's ratio of pay zone 0.2 Dimensionless ν3 Poisson's ratio of bottom layer 0.29 Dimensionless αf,αg Shear failure parameter 0 Dimensionless Af, Ag Shear failure parameter 393.8 Bf, Bg Shear failure parameter 1.869 × 10 psi-1 R Cap parameter 1.559 Dimensionless ψ Triaxial extension to compression strength ratio 1.0 Dimensionless X0 initial value of cap hardening parameter X -200.0 Psi W Permissible maximum inelastic volumetric strain 0.21 Dimensionless D1 Cap hardening parameter 5.0 × 10-5 psi-1 D2 Cap hardening parameter 0 psi-2 σ00 initial value of hardening parameter sigma0 0.6 Psi M Biot's modulus 8.0 × 104 Psi kp/μ Kozeny-Carman fitting coefficient 3.5 × 10-2 ft2/psi/day α Parameter for fluid flow 1.0 Dimensionless C1 Creeping parameter 1.0 × 10-4 Dimensionless z1 Creeping parameter 1.4 × 10-4 Day-1 z2 Creeping parameter 1.0 psi-1 Psi -3 Figure 3—Subsidence in the reservoir at different stages. (a) after 20 days; (b) after 200 days; (c) after 300 days; (d) after 785 days. Figs. 3 presents the subsidence contour of ground surface at four time stages. Because top and bottom layers are very stiff compared to the pay zone, the deformation in the bottom layer is very small, and top 8 SPE-182721-MS layer just drops with the pay zone. Maximum subsidence occurs right above the middle wellbore as shown in Fig. 3. Because of oil extraction, pore pressure diffuses with time and the subsidence is becoming larger and larger. The simulation shows that a serious subsidence of the ground surface is up to 13 feet at 785 days as shown in Fig. 3(d). The subsidence of the central point of top ground surface following time is plotted in Fig. 4. Figure 4—Comparison of the subsidence following time at the central point of the top surface: linear elastic and cap plasticity with creeping models In addition, we also run models with linear elasticity and the Drucker-Prager plasticity without cap. The contour of subsidence obtained with linear elastic and Drucker-Prager model is also shown in Fig. 5 for comparison. Because this problem is dominated by compaction, payzone formation doesn't yield for the case where the Drucker-Prager plasticity model is used. Therefore, it is understandable that the simulation results from the linear elastic and Drucker-Prager models are the same. The subsidence of the central top point increases linearly for a very short time and later the subsidence rate approaches zero as shown in Fig. 4. As shown in Fig. 4, the subsidence predicted from the cap model is as 7 times as the prediction from the linear elastic model or the Drucker-Prager model. As discussed in (Bruno et al. 1992), subsidence for Belridge filed was estimated by 30 ft. with 300-400 ft. soft formulation. We see that the use of the linear elasticity model or the Drucker-Prager plasticity model may seriously underestimate the subsidence of reservoir with soft rocks undergoing compaction. Figure 5—Subsidence of the reservoir with linear elastic/Drucker-Prager model at 785 days Fig. 6 shows that plasticity occurs around 3 wells at earlier time and eventually, the formation of whole payzone yields following the elapse of time. Up to 17% volume plastic strain develops around the middle wellbore at 785 day. From Fig. 6(d), the average volume plastic strain resulting from reservoir compaction SPE-182721-MS 9 is up to 5~10%. For compaction simulations, contours on volume plastic strains on payzone layers are able to provide a good check on if the formation yields on the cap portion of the yield surface. Fig. 7 presents a relatively small creep strain contribution to the subsidence as this problem is dominated by cap plasticity. For longer time, the contribution from creeping effect is expected to substantially increase. Figure 6—Plastic volumetric strain near wellbores at different stages. (a) after 20 days; (b) after 200 days; (c) after 300 days; (d) after 785 days. Figure 7—Creeping volumetric strain in the pay zone at different stages. (a) after 200 days; (b) after 785 days; The path of stress and evolution of yield functions at the right corner of the middle well is drawn in meridional space and shown in Fig. 8. In Fig. 8, there are 6 yield surfaces corresponding to 0, 20, 200, 400, 600 and 780 days, respectively. It can be seen that yield happens at this point a little before 20 days for the first time, and it's already in the compression cap part of the yield surface. This specific point keep staying in the plastic state, and stress is always on the cap portion of the yield surface. This demonstrates that our IGAR computer code is able to model the cap plasticity with large deformation. Figure 8—Yield surfaces at 6 different stages and stress path in meridional coordinate system at the right corner of the middle wellbore. 6 yield surfaces is of 0, 20, 200, 400, 600 and 780 days, respectively. 10 SPE-182721-MS Wellbore Damage To further demonstrate the performance of our proposed computational framework, a wellbore model having a square domain with a wellbore in the center shown in Fig. 9 is solved using IGAR code. Uniform pressure F=80 psi is applied on the top and bottom surface of the model. Left and right surfaces are fixed in normal direction. The pore pressure around the circular well is prescribed as p=10 psi. Solid skeleton parameters are listed in Table 2. Fluid parameters are the same with the reservoir model. Initial stresses are σx0= σz0=-5 psi, σy0=-10 psi. Initial pore pressure is p0=50 psi. Figure 9—Geometry of wellbore model Table 2—MATERIAL PARAMETERS FOR WELLBORE DAMAGE MODEL Parameter Description Value Unit E Young's modulus of top layer 5.0 × 103 psi ν Poisson's ratio of top layer 0.2 Dimensionless αf,αg Shear failure parameter 0.25 Dimensionless Af, Ag Shear failure parameter 0 psi Bf, Bg Shear failure parameter 0 psi-1 R Cap parameter 2.0 Dimensionless ψ Triaxial extension to compression strength ratio 1.0 Dimensionless X0 Initial value of cap hardening parameter X -80.0 psi W Permissible maximum inelastic volumetric strain 0.3 Dimensionless D1 Cap hardening parameter 5.0 × 10-2 psi-1 D2 Cap hardening parameter 0 psi-2 σ00 Initial value of hardening parameter sigma0 6.0 psi Fig. 10 shows original mesh and deformed mesh at 2000 days. From Fig. 10 we can see some serious damage around the wellbore. The wellbore is distorted heavily, reduced by nearly 30% in height. Fig. 11 presents the plastic volumetric strain and mean stress at 2000 days. As shown in Fig. 11(b), mean stress is very large at the left and right corner of the wellbore, and it causes large plastic deformation, which can be seen from Fig. 11(a). Therefore, the plasticity at the right corner is dominated by compaction. Mean stress and plastic strains are small at top and bottom corner. The most serious damage tending to occur at left and right corners. SPE-182721-MS 11 Figure 10—Comparison of mesh around the wellbore. (a) original mesh; (b) mesh at 2000 days. Figure 11—Contour of plastic volumetric strain and mean stress after 2000 days for the wellbore model. (a) plastic volumetric strain; (b) mean stress. Conclusions In this work, we propose a comprehensive computational framework for modeling reservoir compactions. First, to capture serious compaction effect we apply a cap plasticity model rather than traditional DruckerPrager plasticity models. Second, to take into account the long term effect on compaction prediction, a creep model is integrated with cap plasticity models. Third, to achieve fast convergent rates for highly nonlinear problems we have derived a consistent material integrator and incorporated it into a general finite element code for finite deformation problems, for both smooth rate independent cap plasticity and creep models. Finally, a fully coupled geomechanics and fluid flow finite element framework have been adopted. Our comprehensive study on reservoir compaction includes coupled fields, cap plasticity, creep effect, complex payzone geometries with multiple wells, and large deformation. Moreover, Comparisons among 3 rock material models including linear elasticity, Drucker-Prager plasticity, and cap plasticity with creeping behaviors are also performed. This work is expected to serve as a fundamental guide for modeling serious reservoir compactions by incorporating correct physics as well as adopting robust and accurate computational frameworks. References Baade, R., L. Chin and W. Siemers, "Forecasting of Ekofisk reservoir compaction and subsidence by numerical simulation." Journal of Petroleum Technology 41(07): 723 - 728(1988). Biot, M. A., "General Theory of Three-Dimensional Consolidation." Journal of Applied Physics 12(2): 155–164(1941). Bruno, M., "Subsidence-induced well failure." SPE Drilling Engineering 7(02): 148–152(1992). Bruno, M. and C. Bovberg. Reservoir compaction and surface subsidence above the Lost Hills Field, California. The 33th US Symposium on Rock Mechanics (USRMS), American Rock Mechanics Association(1992). 12 SPE-182721-MS Chin, L., R. R. Boade, J. H. Prevost and G. H. Landa, "Numerical simulation of shear-induced compaction in the ekofisk reservoir." International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 30(7): 1193–1200(1993). De Rouffignac, E., P. Bondor, J. Karanikas and S. Hara. Subsidence and well failure in the South Belridge diatomite field. SPE Western Regional Meeting, Bakersfield, California, Society of Petroleum Engineers(1995). Dean, R. H., X. Gai, C. M. Stone and S. E. Minkoff, "A comparison of techniques for coupling porous flow and geomechanics." SPE Journal 11(01): 132–140(2006). Desai, C. S., S. Somasundaram and G. Frantziskonis, "A hierarchical approach for constitutive modelling of geologic materials." International Journal for Numerical and Analytical Methods in Geomechanics 10(3): 225–257(1986). Fossum, A. F. and J. T. Fredrich. Cap plasticity models and compactive and dilatant pre-failure deformation. 4th North American Rock Mechanics Symposium, Seattle, Washington, American Rock Mechanics Association(2000). Foster, C. D., R. A. Regueiro, A. F. Fossum and R. I. Borja, "Implicit numerical integration of a three-invariant, isotropic/ kinematic hardening cap plasticity model for geomaterials." Computer Methods in Applied Mechanics and Engineering 194(50–52): 5109–5138(2005). Fredrich, J., J. Arguello, B. Thorne, W. Wawersik, G. Deitrick, E. De Rouffignac, L. Myer and M. Bruno. Three-dimensional geomechanical simulation of reservoir compaction and implications for well failures in the Belridge Diatomite. SPE Annual Technical Conference and Exhibition, Denver, Colorado, Society of Petroleum Engineers(1996). Gai, X., R. H. Dean, M. F. Wheeler and R. Liu. Coupled geomechanical and reservoir modeling on parallel computers. SPE Reservoir Simulation Symposium, Houston, Texas, Society of Petroleum Engineers(2003). Hansen, K., M. Prats and C. Chan, "Modeling of reservoir compaction and surface subsidence at South Belridge." SPE Production & Facilities 10(03): 134–143(1995). Keszthelyi, D., D. K. Dysthe and B. Jamtveit, "Compaction of North-sea chalk by pore-failure and pressure solution in a producing reservoir." Frontiers in Physics 4(4)(2016). Lewis, R. W. and B. A. Schrefler. The finite element method in the deformation and consolidation of porous media. New York, NY, John Wiley and Sons Inc.(1987). Liu, R., G. Lin and G. Bhashyam, "Transformation of the Sandler and Rubin Nonsmooth Cap Model to the Pelessone Smooth Cap Model." Journal of engineering mechanics 136(5): 680–685(2009). Liu, R., M. F. Wheeler, C. N. Dawson and R. Dean, "Modeling of convection-dominated thermoporomechanics problems using incomplete interior penalty Galerkin method." Computer Methods in Applied Mechanics and Engineering 198(9– 12): 912–919(2009). Liu, R., M. F. Wheeler, C. N. Dawson and R. H. Dean, "On a coupled discontinuous/continuous Galerkin framework and an adaptive penalty scheme for poroelasticity problems." Computer Methods in Applied Mechanics and Engineering 198(41–44): 3499–3510(2009). Motamedi, M. H. and C. D. Foster, "An improved implicit numerical integration of a non-associated, three-invariant cap plasticity model with mixed isotropic–kinematic hardening for geomaterials." International Journal for Numerical and Analytical Methods in Geomechanics 39(17): 1853–1883(2015). Oden, J. T. The Finite Elements of Nonlinear Continua. New York, United States, McGraw-Hill(1972). Pelessone, D. A modified formulation of the cap model. Gulf Atomics Report GA-C19579 to the Defense Nuclear Agency(1989). Rattia, A. J. and S. Ali. Effect of formation compaction on steam injection response. SPE Annual Technical Conference and Exhibition, San Antonio, Texas, Society of Petroleum Engineers(1981). Rinehart, A., S. T. Broome and T. Dewers. Mechanical variability and constitutive behavior at reservoir conditions of the Lower Tuscaloosa Formation supporting the SECARB Phase III CO2 Injection Program at Cranfield Site. Albuquerque, NM, Sandia National Laboratories(2014). Sandler, I. S., F. L. Dimaggio and G. Y. Baladi, "Generalized cap model for geological materials." Journal of Geotechnical and Geoenvironmental Engineering 102(7): 683–699(1976). Schwer, L. E. and Y. D. Murray, "A three-invariant smooth cap model with mixed hardening." International Journal for Numerical and Analytical Methods in Geomechanics 18(10): 657–688(1994). Settari, A. and F. Mourits, "A coupled reservoir and geomechanical simulation system." SPE Journal 3(03): 219–226(1998). Settari, A. and D. A. Walters, "Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction." Spe Journal 6(03): 334–342(2001). Sun, W., Q. Chen and J. T. Ostien, "Modeling the hydro-mechanical responses of strip and circular punch loadings on water-saturated collapsible geomaterials." Acta Geotechnica 9(5): 903–934(2014). Vudovich, A., L. Y. Chin and D. R. Morgan, "Casing deformation in Ekofisk." Journal of Petroleum Technology 41(07): 729 - 734(1989).