Parallel High-Performance Computing in Geomechanics with Inner/Outer Iterative Procedures

Cost-effective improvements in the technology needed to develop and manage reservoirs in challenging environments require an increase in our understanding of geomechanical behavior. The local stresses relevant to reservoir-scale processes are strongly affected by stratigraphic and structural features at that scale, and the relationship between hydrocarbon production and the mechanical behavior of the reservoir and/or overburden can be complicated and difficult to discern from field data directly. Numerical simulation provides a means to achieve critical insight into the behavior of complex geosystems and advance understanding of both the subsurface environment before drilling as well as the relationship between fluid flow and geomechanical behavior during production. This paper reviews our recent work aimed at utilizing large-scale geomechanical simulation as a reservoir management tool. We first describe the constitutive models developed specifically for two important classes of geomaterials that are implemented in the quasi-static largedeformation finite element code JAS3D that we use in our work. We then describe several field cases where we apply nonlinear finite element modeling to key problems in reservoir mechanics. The first field cases involve historical geomechanical simulations of primary and secondary recovery at the Belridge Diatomite and Lost Hills fields located in California's San Joaquin Valley. In this work, we apply nonlinear finite element modeling to investigate the causes of well casing damage experienced by the field operators and to identify mitigation strategies. Next, we describe an application that addresses potential well integrity issues associated with sub-salt and near-salt deepwater Gulf of Mexico reservoirs. In this work, we analyze hole closure behavior and quantify loading on casings for wells that penetrate thick salt formations to ensure that wells are designed to withstand salt loading over a service lifetime of 20-30 years. The field studies illustrate the difficult issues encountered in practical applications of large-scale three-dimensional nonlinear finite element geomechanical modeling and suggest areas where research advances could be beneficial. These areas include: development of solids models and finite element meshes from disparate geologic and/or numerical data, development of realistic constitutive models, robust and efficient implementation of these material models in finite element codes to achieve reasonable solution times, experimental rock mechanics data and the natural heterogeneity of geosystems, implementation of far-field (tectonic) stresses and stress initialization, and integration of geomechanical modeling results with other analysis tools.


INTRODUCTION
The technology needed to develop and manage reservoirs economically in challenging environments requires an increase in our modeling ability such that geomechanical behavior can be predicted accurately.While the far-field in situ stress state results from plate-scale tectonic motions, stratigraphic and structural features produce local stresses that strongly affect reservoir-scale processes.Furthermore, the relationship between hydrocarbon production and the mechanical behavior of the reservoir and/or overburden can be complicated and difficult to discern from field data directly.Posing 'what if' questions through numerical simulation provides a means to achieve critical insight into the behavior of complex geosystems and advance our understanding of both the subsurface environment before drilling, as well as the relationship between fluid flow and geomechanical behavior during production.While twodimensional idealized finite element geomechanical analyses are useful for advancing generic understanding of reservoir mechanics, direct application to specific field situations more typically requires development of large-scale threedimensional models.With recent advances in computational capabilities, as well as concurrent improvements in constitutive modeling of pressure-sensitive geomaterials, it is now possible to apply large-scale geomechanical finite element modeling to a variety of key problems in reservoir mechanics.
This paper reviews our recent work aimed at utilizing geomechanical simulation as a reservoir management tool.We first describe briefly the three-dimensional quasi-static nonlinear finite element code JAS3D that we apply in this work and the constitutive models developed specifically for two important classes of geomaterials.We then describe several field cases where we apply nonlinear finite element modeling.The first field cases involve historical geomechanical simulations of the Belridge Diatomite and Lost Hills fields located in California.In this work, nonlinear finite element modeling is applied to investigate the causes of well casing damage experienced by the field operators and to identify mitigation strategies.The approach includes the development of large-scale three-dimensional geomechanical models (hundreds of thousands of finite elements) and unidirectional coupling between three-dimensional black oil reservoir models and geomechanical models.Next, we describe a field application that addresses well integrity issues associated with sub-salt and near-salt deepwater Gulf of Mexico (GoM) reservoirs.In this work, we analyze hole closure behavior caused by salt creep and quantify loading on casings for wells that penetrate thick salt formations to ensure that wells are designed to withstand salt loading over a service lifetime of 20-30 years.Besides directly demonstrating both the feasibility and benefits of large-scale three-dimensional non-linear finite element geomechanical modeling, the field cases also illustrate the difficult issues encountered in practical applications.We close by discussing these issues and the approaches that we have taken to circumvent some of these difficulties.We also suggest areas where research advances would be beneficial.

JAS3D
The finite element code JAS3D (Blanford, 1998), based on massively parallel distributed memory multiprocessor computers, is ideally suited for large-scale geomechanics analyses.The code utilizes explicit quasi-static solvers developed specifically under US Department of Energy (DOE) Defense Programs to handle very large nuclear weapons problems with complex material constitutive models and contact surfaces.The use of parallelized iterative solvers and the extensive experience gained with the use of nonlinear material models provide the base technology that offers the capability for efficient solution of large, geometrically complex problems.Commercially available finite element technology generally cannot handle models of the size that are required to capture adequately the complex geology that is associated with hydrocarbon reservoirs in specific field settings of interest.
JAS3D is capable of predicting the mechanical response of soil and rock over the spatial and temporal scales of interest to the oil and gas industry, including both the geologic scale (hundreds to thousands of kilometers and millions of years) as well as the reservoir scale (tens of kilometers and tens of years).It was envisioned that such a capability could be a valuable tool in reducing risks on prospects, as well as for planning initial development and subsequent management of the reservoir over the course of the field's production lifetime.Specific examples of these applications include: optimizing well trajectories by consideration of the subsurface stress field, predicting subsidence caused by pore pressure drawdown and reservoir compaction, designing (and not over-designing) well casings, and predicting the evolution of reservoir stresses during production.Several applications demonstrate that JAS3D can handle the mechanics applications that are critically important to the exploration and production activities of the oil and gas industry (Fredrich et al., 1996(Fredrich et al., , 2000(Fredrich et al., , 2001;;Minkoff et al., 1999;Willson et al., 2002).

Constitutive Models for Geomaterials
The constitutive behavior of rocks differs fundamentally from the behavior of engineering materials such as metals.Of primary importance for reservoir mechanics is that geomaterials such as sandstone and chalk are pressure sensitive; that is, their strength and deformation vary as a function of mean stress.This has great implications for developing an understanding of the geomechanical response of a hydrocarbon reservoir to production with accompanying changes in pore pressure.
While insignificant from the point of view that salt bodies are generally fully dense in the Earth and therefore have negligible storativity, the constitutive behavior of salt is nevertheless of considerable significance to reservoir mechanics.Besides being central to the geologic evolution of important oil and gas provinces such as the Gulf of Mexico, salt bodies affect the local stress state, and also present challenges during drilling through salt that overlies hydrocarbon reservoirs located at greater depths.Polycrystalline salt exhibits rather unique features in that its deformation behavior is similar to the deformation behavior of rocks in the brittle regime and of metals in the ductile regime.
Through work performed for the US DOE, Sandia National Laboratories has developed a library of constitutive models for geomaterials.In particular, a constitutive model for salt was developed to model the complex behavior of natural salt occurring at the Waste Isolation Pilot Plant (WIPP) site, the US DOE's repository for transuranic waste generated by US national defense programs.Sandia's expertise in modeling salt was then extended further through its work as scientific advisor to the US DOE for the Strategic Petroleum Reserve, a series of underground storage sites located in salt domes along the US Gulf Coast.In the mid-1990s, an effort was initiated to expand the non-salt geomaterials library beyond conventional Drucker-Prager and Mohr Coulomb models to include development of a class of material models specifically applicable to porous geomaterials.Besides being necessary for modeling reservoir behavior, this development was simultaneously motivated by the needs in defense programs for modeling the penetration of projectiles into hard materials such as rock and concrete and assessing our ability to perform vulnerability and survivability assessments of hardened deeply buried targets.There were also needs to model powder compaction used in the manufacture of modern sophisticated weapons components.The two main classes of constitutive geomaterial models mentioned here are described below in more detail.

Compacting and Dilating Geomaterials
Two types of behavior may be identified on the basis of rock porosity for the range of mean stresses encountered in geosystems of engineering interest (e.g. oil and gas reservoirs, nuclear waste repositories, buried targets, and depleted reservoirs for possible use for subsurface sequestration of greenhouse gases).Low porosity rocks exhibit shear localization.High porosity rocks exhibit more varied deformation, including dilation, shear localization, compaction, and compaction localization, depending on the mean stress.At low mean stresses, porous geomaterials fail by shear localization, and at higher mean stresses, they undergo strainhardening behavior.Cap plasticity models (e.g., DiMaggio and Sandler, 1971) attempt to model this behavior using a pressure-dependent shear yield and/or shear limit-state envelope with a hardening or hardening/softening elliptical end cap to define pore collapse.
While the traditional cap plasticity models describe compactive yield and ultimate shear failure, difficulties arise when the rock behavior involves a transition from compactive to dilatant deformation that occurs well before the shear failure or limit-state shear stress is reached.Laboratory rock mechanics tests on porous sandstones and carbonates reveal precisely this behavior.To circumvent this limitation, a continuous surface plasticity model was developed (Fossum and Fredrich, 2000a) to predict compactive and dilatant prefailure deformation.During loading the stress point can pass freely through the critical state point separating compactive from dilatant deformation so that the predicted volumetric strain goes from compactive to dilatant without the use of a non-associated flow rule.

Traditional Cap Plasticity Model
The constitutive equations referred to as cap plasticity models have been in existence for at least 45 years (Drucker et al., 1957).Traditionally, these models comprise two independent continuously differentiable yield functions (Fig. 1a).The inner envelope of the two surfaces corresponding to the two smooth yield functions generates the yield surface that defines the boundary of the elastic region.Corners appear at the intersections of the two yield surfaces.In the case of flow, some workers accept that the plastic strain increment is not specified uniquely but is only required to lie between the normals to the yield surface at adjacent points on the yield surface.Other workers attempt to eliminate this indeterminacy by introducing a third yield function that serves to connect the first two yield surfaces continuously.In the case of hardening behavior, the plastic strain increment comprises the sum of the plastic strain increments determined independently from each of the yield functions.The inner envelope of the surfaces must be convex requiring that each individual surface be convex.The cap surface is required to intersect the shear-failure surface at the point of horizontal tangency to the cap ellipse.Thus, the shear failure surface represents the critical state line separating dilatant from compactive inelastic volumetric deformation.The model 426 has a number of shortcomings.As mentioned above, there is an indeterminacy of flow direction at the intersection of the shear-failure surface and the cap-hardening surface.From a numerical standpoint, it is often found that an excessive amount of solution time is spent iterating to find the intersection of the shear failure surface with the cap hardening surface.From a physical viewpoint, the requirement that the cap hardening surface intersect the shear failure surface at the point of horizontal tangency to the cap ellipse precludes pre-failure dilatant deformation, contradicting experimental observations.Also, a three-invariant model is often required for rocks and similar brittle materials.

Continuous Surface Plasticity Model
To circumvent the shortcomings of the model described in the previous section, a variation of a technique described by Schwer and Murray (1994) is used to smooth the intersection of the failure and cap hardening surfaces using a Pelessone function (Pelessone, 1989).The resulting continuous surface yield function (Fig. 1b) incorporates shear kinematic hardening, isotropic cap hardening, Lode angle dependence of yield, and is sufficiently general to predict pre-failure dilatant deformation.The constitutive equations are described in detail by Fossum and Fredrich (2000a).More recently, the model has been expanded to include anisotropy (more specifically, transverse isotropy) in both the elastic and inelastic regimes, and strain-rate sensitivity for high strain rate applications.Measured and predicted volumetric strains as a function of axial stress for an unconfined compression test performed on diatomaceous sand from the Lost Hills oil field, California, for the traditional cap plasticity model.The traditional cap plasticity model cannot predict the transition from compaction to dilation that occurs well before peak stress is reached (1 psi = 6894.8Pa).
Two features are of note: all subsequent loading surfaces and the initial yield surface are convex, and the plastic strain increment vector is normal to the loading surface.These conditions are consistent with Drucker's hypothesis for a work-hardening material, namely that useful net energy over and above the elastic energy cannot be extracted from the material, and the material is said to be stable.The advantages of the continuous surface plasticity model are illustrated by comparing the predicted behavior for the traditional (Fig. 2) versus new (Fig. 3) model as applied to several data sets for a single reservoir formation.The data are from conventional triaxial compression tests, and the material is a sand unit (D) of the Upper Etchegoin Formation (late Pliocene to early Miocene).The material is a very high porosity unconsolidated diatomaceous sand with a bulk density of ~ 1750 kg/m 3 (1.75 g/cc) and is one of the two reservoir formations at the Lost Hills oil field located in the San Joaquin basin of central California.
Implementation of the material model in the finite element code JAS3D is described in Fossum and Fredrich (2000a).To benchmark the performance of the new model, a largescale three-dimensional numerical simulation was performed for direct comparison with a previous simulation described by Fredrich et al. (2000) that used the traditional model.The computational time using the newly implemented model is one-quarter that required by the traditional model.As discussed above, this dramatic improvement is a direct consequence of the use of a continuous yield surface.

Implications of Compaction Behavior for Reservoir Mechanics
As a reservoir is produced, the pore pressure typically depletes, leading to an increase in effective stress.For weak sandstones or chalks, the increase in effective stress can be sufficient that the stress state loads the compactiondominated region of the yield surface.As discussed by Segall (1989), this inelastic deformation in the reservoir that results from pore pressure depletion affects the stress state in the surrounding sideburden and overburden (Fig. 4).Contraction of the reservoir rock in the vertical direction is accommodated partially by subsidence of the free surface, whereas contraction in the horizontal direction is resisted by the surrounding rock that is pulled towards the reservoir.The resulting stress state can be highly perturbed as well as locally variable.The vertical strains may be sufficient to induce compressional casing failures, and in the presence of pre-existing faults or weak bedding planes, the local changes in stress may also be sufficient to cause shear, or dog-leg, casing failures.Well casing damage induced by formation compaction has occurred in reservoirs in the North Sea, the Gulf of Mexico, California, South America, and Asia.As discussed below, we apply nonlinear finite element modeling to examine this issue for two compacting reservoirs located in California's San Joaquin Valley.

Viscous Rock Salt
Polycrystalline salt exhibits deformation behavior that is similar to the deformation behavior of rocks in the brittle regime and of metals in the ductile regime.Salt creeps under any deviatoric stress and the creep rate is a strong function of the deviatoric stress.For a given mean stress, it takes less energy to deform salt in pure shear than it does in either triaxial compression or triaxial extension.For mean normal stresses typically less than about 5 MPa (725 psi) it will dilate with time upon application of a deviatoric stress, but it will flow with constant volume above a mean stress of ~5 MPa (725 psi).The brittle failure of salt, defined as the shear stress at peak load, is pressure sensitive and a function of plastic strain.
The creep of salt involves either two or three creep regimes.For confining pressures typically less than 5 MPa (725 psi), salt specimens subjected to a constant stress state will pass through three stages.In the first stage, called primary creep, the strain rate begins with a very high rate and decreases to a constant rate.The time during which the specimen deforms at constant rate is called secondary creep.The third stage is called tertiary creep in which the strain rate increases until failure occurs.For confining pressures typically above 5 MPa (725 psi) only the primary and secondary stages are evident.
The constitutive model developed at Sandia for WIPP and SPR analyses is called the Multimechanism Deformation Coupled Fracture (MDCF) model (Munson and Dawson, 1979;Munson et al., 1990;Chan et al., 1994).This model is formulated by considering individual mechanisms that include dislocation creep, shear damage, tensile damage, and damage healing.The MCDF model, which is implemented in Sandia's finite element codes, is significant because it is based, as much as possible, on first principles (i.e., identified   deformation mechanisms), as opposed to a phenomenological model that is based on macroscopic observations.This aspect is important because the former approach is required to extrapolate beyond the conditions used in the laboratory to determine the constitutive parameters for the model.While a phenomenological model may fit experimental data very well, it is not necessarily a good predictor for conditions of pressure, temperature, and time outside of the laboratory experiments.WIPP salt, although from a bedded formation, can be considered representative of the known deformation behaviors of US Gulf Coast salt domes as illustrated in Figure 5. Figure 6 reveals some of the following important characteristics that distinguish the constitutive behavior of salt from the porous geomaterials described above: -the plastic incompressible (isochoric) nature of the deformation at mean stresses greater than 5 MPa (725 psi); -the independence of creep rate on mean stress, but high dependence on maximum shear stress; -the high temperature dependence of the creep rate; -that salt tends to heal under elevated mean stress.Based on validation studies at the WIPP site and on observations of SPR cavern behavior compared with model studies, it is concluded that the current state-of-the-art in salt modeling is adequate to model the deformations and stress fields of underground structural configurations for the ranges of temperature and stresses of interest to WIPP and SPR.Although the conditions expected for deepwater GoM field development are more severe than those encountered at the WIPP or SPR sites, they nonetheless fall within the ranges of stress and temperature for which the mechanistically based salt models have been developed.Willson et al. (2002) present laboratory creep tests performed on rotary sidewall cores taken from a deepwater GoM salt body adjacent to the Mad Dog field.The laboratory tests reveal this deepwater salt to be a relatively "strong" salt (see Fig. 5), and, moreover, show the creep response of the Mad Dog salt to be nearly identical to the creep behavior of the Bayou Choctaw salt (Fig. 7).Besides representing the deepest salt cores ever subjected to controlled laboratory creep tests, this work is also significant because it suggests that the known behavior of SPR Gulf Coast salts (Munson, 1999) can be expected to typify the behavior of deepwater GoM salts.

Implications of Salt Constitutive Behavior for Reservoir Mechanics
Figure 8 shows steady-state strain rates calculated for WIPP salt using the MCDF constitutive model for the range of shear stresses and temperatures expected in deepwater GoM field developments.Also shown are stress profiles as a function of depth from mudline (nominal 4000 ft water depth, 1 ft = 0.3048 m) that are expected to bound the conditions encountered in the deepwater GoM.A generic overburden profile was used to calculate the vertical stress profile for normally pressured formations, and two empirical relationships determined from field data were used to compute the magnitude of the minimum horizontal stress (SM Willson, Private Communication, 2001).For both profiles, the stress difference used for the calculations equals the overburden stress minus the minimum horizontal stress, and the in situ temperatures follow a relationship of 38°F at mudline, with a temperature gradient of 1.45°F per 100 ft depth.
It is critical to recognize that the strain rates shown in Figure 8 are creep rates.The significance of this is that creep rates correspond to fixed or constant stresses.Since stress is a function only of the elastic strain (i.e., not the inelastic or total strain), this implies that the elastic strain rates are zero.However, this is not the case in the subsurface where the stresses in the region of interest are affected by the constraint imposed by the overburden, sideburden, and underburden formations.
In a laboratory creep test, a specimen is loaded to a certain stress level and then the stresses are held constant while the specimen is allowed to deform with time.Except for the elastic load-up strain, all subsequent strain is inelastic, and since the elastic strain is zero, this also represents the total strain.Thus, the creep test represents a structural problem that is "statically determinate", meaning that the stresses can be determined from equilibrium equations alone.However, this is generally not the case in the subsurface; instead the full set of equilibrium equations, compatibility equations, and constitutive equations must be solved with the necessary initial and boundary conditions.This is a much more Comparison of steady-state creep behavior of deepwater GoM salt and Bayou Choctaw SPR salt (after Willson et al., 2002) (1 psi = 6894.8Pa).complicated problem than the simple laboratory test, and accurate prediction of the deformation rates in a real field setting requires application of numerical methods such as the finite element method.Only with numerical methods can one simultaneously solve the equilibrium equations, compatibility equations, and constitutive equations subject to the imposed initial and boundary conditions.
If a salt body in the subsurface is suddenly subjected to a significant stress difference, this stress difference will produce instantaneously an inelastic strain rate.Without transient effects, this inelastic strain rate will equal the strain rate corresponding to the stress difference as shown in Figure 8.However, in the subsurface the salt will be constrained and the elastic strain rate will not be zero.The total strain rate will be the sum of the elastic strain rate and the inelastic strain rate.Because of constraint, the total strain rate may be small or even close to zero.If the inelastic strain rate were high but the total strain rate were constrained to be small, then the elastic strain rate would have the opposite sign of the inelastic strain rate such that the sum of the elastic and inelastic strain rates equals the total strain rate.This process results in a stress change.So, while very high stress differences in the subsurface may produce very high inelastic strain rates in the salt, the in situ stresses will likely change very rapidly as well.The only way to determine what the total strain rate is in the salt, and thus how much the salt moves, is to solve the complete set of equilibrium, compatibility, and constitutive equations with the appropriate initial and boundary conditions.
Because of the requirement for a salt body to be in equilibrium and to maintain continuity with the surrounding formations, the stress state at the interface with the salt body may be highly complex and variable.Consider the simple case of a spherical salt body that is surrounded by a non-salt formation with different principal vertical and horizontal stresses.The stresses in the salt will relax and ultimately reach a hydrostatic stress state, as long as the far field stresses are not perturbed, and as long as the surrounding formation is constrained.But because of the stress relaxation in the salt body, the stress state in the non-salt formation immediately surrounding the salt will be different from the far field stress state to maintain equilibrium and continuity.Generally, the principal stresses will rotate, or equivalently, shear stresses will appear in the horizontal and vertical planes.This is illustrated for the simple case of a spherical salt body in Figure 9 that shows the von Mises stress, defined as the square root of one-half of the sum of the squares of the principal stress differences.These stress rotations have implications for the long-term integrity of well casings at the interface between salt and non-salt formations.Further in the paper, we discuss another salt-related application that considers the effect of salt creep on the loads experienced by a through-salt well casing over the anticipated well service lifetime.10 -15 s -1 10 -14 s -1 10 -13 s -1 10 -12 s -1 10 -11 s -1 10 -10 s -1 10 -9 s -1 10 -8 s -1 10 -7 s -1 10 -6 s -1 10 -5 s -1 10 -4 s -1 10 -3 s -1 10 -2 s -1 10 -1 s -1 10 -0 s -1 10 +1 s -1 10 +2 Contour map of steady-state creep strain rates for salt, in the absence of geometrical and structural constraints.The two lines bound the expected stress and temperature conditions for deepwater GoM field developments (see text).

Reservoir Compaction and Casing Damage at the Belridge Diatomite and Lost Hills Oil Fields, California
The shallow diatomaceous oil reservoirs located in Kern County, California, are susceptible to depletion-induced compaction because of the high porosity (45-60%) and large vertical extent of the producing formation.The Belridge Diatomite and Lost Hills fields are two of several significant fields located in the San Joaquin basin.Other major fields include Cymric, Elk Hills, Buena Vista, and Midway-Sunset.
The producing structures are a northwest-southeast trending doubling-plunging anticline that roughly parallel the San Andreas fault located to the west.More than a thousand wells at the Belridge and Lost Hills fields have experienced severe casing damage as the reservoir pressures have been drawn down.Fredrich et al. (1996Fredrich et al. ( , 2001) ) present casing damage statistics for the Belridge and Lost Hills fields, respectively.Although there are faults in the field area, these faults are not known to be active (Fredrich et al., 2000).Moreover, the casing damage statistics (Fig. 10) as well as three-dimensional visualizations of casing damage (Myer et al. 1996) show that damage is focused along bedding planes located at the top of the diatomite reservoir and in the overburden, rather than along shear-oriented faults (Fredrich et al., 2000(Fredrich et al., , 2001)).
Analyses of field-wide databases indicated that plane-strain approximations could not capture the locally complex production, injection and subsidence patterns, and motivated development of three-dimensional geomechanical models to evaluate potential well failure mechanisms.Historical simulations were performed for two independently-developed areas (Sections 33 and 29) at the Belridge Diatomite field, and for the core area at Lost Hills.
In each of the three applications, time-dependent reservoir pressures are derived from three-dimensional black oil reservoir simulations and applied as loads in three-dimensional nonlinear finite element simulations.The goal of the modeling effort is to examine the role of local gradients in effective stress (that result from gradients in pore pressure induced by aggressive production and injection) in affecting well casing integrity.The close well spacing at these fields (less than a hundred feet in some cases) and resultant high degree of well interaction necessitate the development of very large-scale models.The geomechanical models covered spatial areas on the order of thousands of feet and included hundreds of wells.The models (Fig. 11) are meshed from structural data corresponding to marker horizons within the reservoir and stratigraphic boundaries within the overburden.Contact surfaces are defined along bedding planes at multiple locations in the overburden and at the interface with the reservoir, as inferred by the field data that suggest that these features localize casing shear (Fig. 10).Geomechanical units are defined based on combined analysis of rock Contours of von Mises stress around a spherical salt body that is constrained in a uniform overburden with unequal vertical and horizontal stresses.The stress differences in the salt body approach zero because of the inability of the salt to maintain deviatoric stresses.Consequently, the stress state in the near-field takes up the load that the salt sheds, i.e., the near-field stresses are perturbed from the far-field conditions as required to satisfy the equilibrium conditions.The stress difference in the far-field is as defined by the initial in situ stress state applied.mechanics test data for different intervals within the producing formation, geologic marker data, and lithologic correlation sections.
The structural data are extrapolated outside of the reservoir area to construct two flank regions on either side of the produced region (each with an area equal to the produced area).For the relatively shallow reservoir depths at Belridge and Lost Hills, it was deemed adequate to extend the lateral model boundaries by this factor of × 2 in either direction.(For each of the three geomechanical models, the negligible displacements in the outer parts of the flank regions suggests that the boundaries are indeed sufficiently removed from the reservoir to eliminate boundary effects, e.g., see Figure 6 of Fredrich et al., 2001).In addition, contact surfaces are defined at the top of the reservoir and within the overburden based on the casing damage statistics.
To assess the predictive capability of the model, the predicted surface subsidence is compared with observed historical surface deformations derived from survey monument data at several locations across the model area with generally good agreement (see Fredrich et al., 2000Fredrich et al., , 2001)).The two geomechanical models developed for Belridge illustrate the dramatic consequences of different reservoir management practices on long-term well integrity.Histogram showing the depth of well casing damage at the Belridge Diatomite field, 1984-1994.Note that to correct for the variation in structural depth across the field, depth of casing damage is scaled so that the interface between the top of the reservoir and overburden is equal to zero, and the production zone is shown with shading.The field data show that casing damage is localized on two prominent bedding planes that are located at the interface between the reservoir and overburden, and at a second location some 300 ft higher in the overburden.
The Section 33 geomechanical model included 197 wells and covered 18 years of production history; field data indicated that nearly a hundred of these wells had experienced major casing damage.The Section 29 model included 85 production and injection wells and comprised 19 years of production history; however, only about a dozen of the wells experienced casing damage.Whereas the Section 33 simulation predicts that large regions undergo incremental shear displacements of 1/2 ft/year over repeated successive years, predicted shear deformations for the Section 29 simulation reach 1/2 ft/year only in a very small number of instances with extremely restricted lateral extent (e.g., see Fig. 6 of Fredrich et al., 2000).Comparison of the results for Sections 33 and 29 thus provides a direct global measurement of their relative proclivity for well casing damage and corroborates the vastly different well failure histories observed for the two sections.
An important result of the Section 33 model relates to the lack of uniformity in the predicted well deformations that show remarkable individuality in their response to production and injection.This reflects the highly localized nature of the spatially and temporally varying pore pressure field (and thus effective stress) that surrounds wells, the influence of nearby producers and injectors on well response, and the central role of the local stress changes in driving casing damage.The horizontal and vertical displacement fields at one example well location are shown in Figure 12 and reveal several important features.First, the large change in mechanical properties at the interface between the reservoir and overburden formations coupled with the sliding contact surface introduces a major discontinuity in the horizontal displacement field with depth.Second, the predicted horizontal x-displacements are small until 10 years, which is just before the onset of the waterflood and conversion of this well from production to injection.This well's conversion from production to injection coincides exactly with the development of severe shearing displacements at the well location.This shows clearly the waterflood that was implemented in part to reduce the well failure risk in some cases amplified the horizontal shearing displacements that were experienced by wells, because the water injection in effect exacerbates the pore pressure gradients induced by production.Displacement in the horizontal y-direction (middle plot) is substantially larger than that in the xdirection, and this largely reflects the formation of a subsidence bowl as material flows inwards in response to pressure depletion.The right plot shows that the maximum vertical displacement occurs at the top of the diatomite reservoir, rather than at the surface.Thus, differential compaction caused by production from the underlying reservoir puts the overburden in a state of extension in comparison to the pre-production condition.
The differential compaction of the reservoir formations causes rotations in the relative directions of the principal stresses.As production progresses, the minimum principal stress, that is initially oriented in the horizontal plane, rotates to the vertical direction in the upper intervals of the reservoir (Fig. 13).The rotations have important implications for the orientation of hydraulic fractures during in-fill drilling, and the predicted stress rotations are corroborated by field data that indicate the occurrence of horizontally oriented hydraulic fractures during infill drilling (Wright and Conant, 1995).
The geomechanical modeling at Belridge and Lost Hills demonstrates clearly both the utility and practicality of applying three dimensional geomechanical modeling to understand and predict behavior at both the well and reservoir scale during production.Ultimately, it is highly desirable to apply geomechanical simulation to assess quantitatively the effects of different producer to injector ratios, well density (spacing), the relative timing of infill wells, and initiation or expansion of secondary recovery methods such as waterflooding on well failure potential and subsidence.

Well Casing Integrity Analyses for Sub-Salt and Near-Salt Deepwater GoM Fields
The deepwater GoM is the most active deepwater region in the world and currently provides some of the greatest challenges in scope and opportunity for the industry.The deepwater GoM is estimated to contain undiscovered recoverable resources of at least ~13 Bbbl of oil equivalent (boe), and is known to harbor some exceptional reservoirs such as the recent Crazy Horse discovery (recently renamed Thunder Horse) in over 6000 ft of water with estimated recoverable reserves of at least 1 Bboe.However, the complex salt tectonics and extreme water and reservoir depths result in very high development costs, in addition to requiring innovative technology to bring these fields on stream.Integral to the successful economic development of these fields, where the cost of a single well is tens of millions of dollars, is a well lifetime of 20 to 30 years.A majority of the wells for the most significant developments will penetrate potentially considerable thicknesses of salt (1000 ft to 6000 ft not being uncommon).Assuring the integrity of these wells over the life of the field is a major drilling engineering challenge, as the consequences of well failures may result in billions of dollars remedial costs and lost production (Willson et al., 2002).While the constitutive behavior of salt is well known as discussed above, our knowledge of the influence of salt deformation at both the well-scale and reservoir-scale is poor.Of primary concern is that the lateral movement of salt may jeopardize the integrity of well casings drilled through salt as a consequence of nonsymmetric loading and induced shears at the bounding formation interfaces.
To address this issue, a major finite element modeling effort was initiated to assess the loading on a well casing over the well service lifetime.Nonlinear finite element analyses are being conducted at both the reservoir and wellbore scales, although here we discuss only the wellbore-scale modeling effort.Included in the analyses is transient creep behavior during an excavation (drilling) phase of approximately oneweek duration, in which the drilling mud weight pressure is imposed on the inner surface of the borehole.To mimic barite settling in the mud in the annulus between the salt and the casing, the mud pressure is then replaced by the appropriate sea water pressure after the one-week drilling phase.Subsequent closure behavior for the internally pressurized casing caused by both transient and steady-state creep are then analyzed over a 20-30 year service lifetime.Initial analyses have considered several casing sizes and stress and temperature conditions appropriate to deepwater GoM reservoirs.
Analyses have been performed to assess salt loading for a cemented versus uncemented annulus, and to evaluate potentially nonuniform loading as a result of imperfect hole 435 Figure 13 Contour plot showing the direction of the minimum principal stress in the uppermost reservoir layer after (left) 10 years and (right) 15 years of production.Initially, the minimum principal stress is in the horizontal plane (red), but with time, it re-orients to the vertical direction (blue).quality (i.e., perfectly circular versus slightly elliptical boreholes).An example of the results for 20 inch casing set in a 24 inch ovalized hole are shown in Figure 14.For these analyses, the stresses in the salt formation are first prescribed as initial stresses, with the vertical and horizontal stresses assumed to be equal, giving an initial lithostatic stress state.Normal displacements on the outer edges of the model, that is extended out approximately fifty times the borehole radius, are prescribed to be zero.The model assumes plane strain along the direction of the well casing, and the well casing is assumed to be perfectly plastic.The analyses reveal several important results.First, for a perfectly circular borehole, loading over the well service lifetime is uniform and not sufficient to induce casing yield for the casing sizes considered in the analyses to date.Ultimately, the magnitude of the loading on the outer casing is equal to the far field lithostatic stress.
Imperfect hole quality greatly impacts the induced loads on the well casing and causes yielding of the casing for stress and temperature conditions expected in some field settings (Fig. 14d, 14e, 14f).Under these circumstances, for the case of an uncemented borehole, plastic hinges develop and the casing eventually ovalizes to the shape of the original noncircular borehole, after which the subsequent deformation is stabilized because the loading becomes uniform with time.The presence of cement serves to distribute the loading more uniformly, resulting in a more stable configuration (Fig. 14g, 14h).The results underscore the critical need for such analyses to design casing strings properly for these field developments.Additional analyses for a specific deepwater GoM field currently under development are described in Willson et al. (2002), which also compares the Sandia JAS3D analyses with results of an analytical model used in the industry.Willson et al. (2002) demonstrates that application of the analytic model greatly overestimates the expected loading on well casings that penetrate salt, and discusses the role of the finite element analyses in minimizing the over-design of through-salt well casings.A more systematic and complete set of analyses is currently being performed with the goal of enabling development of a screening program whereby operators may assess the approximate loading on casings, and based on that, determine the need, or lack thereof, for exact analyses for the specific conditions of interest.

THE CHALLENGES OF REAL FIELD APPLICATIONS
The above examples illustrate how the current state of the art of finite element codes, including formulation of relevant constitutive models, and available computational power, directly enables the application of these tools to the complex problems typical of the oil and gas industry.Nevertheless, several outstanding issues encumber the application to realworld problems and therefore suggest directions for future developments.

Development of Solids Model and Finite Element Mesh
In addition to developing sophisticated finite element software to meet the demands of the nation's defense programs, Sandia has developed the attendant solids modeling tools and three-dimensional finite element mesh generators.These tools are ideally suited to create the complex threedimensional meshes that are typical of engineered structures.
CUBIT is a two-and three-dimensional finite element mesh generation toolkit that is being developed to pursue the goal of robust and unattended mesh generation (Mitchell et al., 2000).CUBIT generates surface and volume meshes for solid model-based geometries that are used for finite element analyses.A combination of techniques including paving, mapping, sub-mapping sweeping and various other algorithms under development are available for discretizing the geometry into a finite element mesh.While CUBIT is specifically designed to reduce the time required to generate all-hexahedral meshes, it also provides the capability to generate hex-dominant and tetrahedral meshes.Geometry creation is accomplished using the geometric primitives and boolean operations or by reading models from a variety of file formats including the ACIS SAT or equivalent file format.
Rarely, however, are oil company partners currently able to provide data in this ideal format.Instead, it is more common that relatively complex scripts must be written to port and/or interpolate the geometry from other applications that are typically provided in ASCII data format.A more significant difficulty is presented by the geometry of geologic structures that tend to be characterized by irregularly shaped layers, including pinch outs.These geometries are not well suited to the methods described above based on Boolean operators and geometric primitives.These aspects can therefore greatly increase the complexity and manpower effort required to build a reservoir-scale finite element mesh.There is considerable need for more flexible methods to generate the irregular geometries typical of geologic systems.

Development of Realistic Constitutive Models
The constitutive models used in geomechanical modeling should be detailed enough to describe the expected ranges of material behaviors without being overly complex.For example, in most cases it is not necessary to describe the underburden materials below the reservoir with anything other than a simple elastic model.Likewise, for unusual classes of geomaterials that are relatively homogeneous such as salt, it may be possible to adequately bound the possible behavior using the existing Sandia databases developed for WIPP and SPR salts.Further, for some problems such as determining the local stress field around salt bodies, analyses only need to consider steady state creep.However other analyses, such as predicting salt loading on a well casing over the field lifetime, require the inclusion of transient creep.For weakly consolidated reservoirs, it is critical to consider the ability of the material model to represent compaction accurately, as well as dilation and shear failure.
Notwithstanding recent developments in constitutive modeling of compacting and/or dilating geomaterials, our fundamental understanding of the constitutive behavior of geomaterials is still not complete.From a fundamental perspective, cap hardening behavior implemented in cap plasticity models is generally physically guided, but not physically based.At present, functional relationships are developed on a phenomenological basis to match available experimental data.It would be desirable to have a model that is micromechanically based, so that it could be truly predictable, but such a task should be tempered by the accuracy called for by the specific application.More significantly, complicating aspects of the constitutive response are often poorly understood from a physical basis, and therefore ignored in general constitutive models.Incomplete knowledge of the material response can have catastrophic consequences, such as experienced at Ekofisk, where the chalk weakening effect was not predicted in the original forward modeling of reservoir compaction following implementation of the waterflood (Nagel, 2002).Another example is the diatomite fields of California, where the wellknown fact that diatomite creeps under static load conditions is not considered in the constitutive models developed to date.Similar behavior is likewise observed for weakly consolidated sands.While it may be appropriate to neglect this time-dependent behavior for semi-quantitative analyses, such effects need to be incorporated if it is truly desired to use geomechanical analysis to guide reservoir development and management decisions quantitatively.

Robust Implementation of Constitutive Models in Finite Element Codes
Ideally, the material model needs to capture the physical details of the material response accurately.However, there is a clear trade-off between the sophistication of the material model versus the additional computational time that will be spent solving the rate forms of the constitutive equations, which can negatively impact the ability to perform multiple simulations for large models.Thus, the practical implementation of the model in a numerical finite element code needs to be considered during the conceptual model development, as well as efficient techniques for its solution.
The complex material models required in geomechanics applications can be represented by an initial value problem associated with the solution of a set of ordinary differential equations of first order.These equations must be integrated accurately in the context of a finite element program.For the integrator to be robust, it must maintain some adaptive control over its own progress thereby making frequent changes to the size of its time step as it advances along the path of its solution.The objective of automatic time-step control is to achieve a stable solution of predetermined accuracy with minimal computational effort.The strain-rate equations for salt creep are particularly difficult to integrate because they are said to be numerically stiff, and special techniques have been developed that result in a robust integrator.Cap plasticity models for compacting and/or dilating geomaterials are significantly more computationally intensive than the simpler Drucker-Prager or Mohr Coulomb material models.However, formulation of a continuous surface plasticity model, as opposed to the two-yield surface model originally developed, reduced the run time for a largescale simulation by 75%.In JAS3D, the continuous surface model is integrated by a forward differencing procedure using the complete flow rule.The multiple yield surface model, on the other hand, requires iteration to locate the intersection of the shear yield surface with the cap surface.

Parameter Estimation, the Natural Heterogeneity of Geologic Formations, and Uncertainty
The constitutive models contain material parameters that naturally must be determined from laboratory rock mechanics data.Of course, it is important that the proper experimental data are available so that the material parameters are well-determined and unique.Fossum andFredrich (1998, 2000b) describe material parameter estimation procedures for cap plasticity constitutive models from experimental rock mechanics test data.The larger issue underlying material parameter estimation is availability of core, and a willingness to invest in a laboratory test program.Decisions regarding how many intervals of the reservoir should be tested for determination of multiple material models within the producing interval, and the number of cores that should be acquired from different regions within the reservoir, should be based upon examination of available geologic and well log data.Likewise, the necessity for including other mechanical features such as faults or joints needs to be based on geologic observations.This in turn may require associating a friction coefficient for the fault, or joint set stiffness properties.
It is likewise critical for coupled reservoir flowgeomechanical modeling that the reservoir flow models adequately represent the fluid pressure distribution and flow field in the reservoir.Besides having again the common factor of material heterogeneity in physical properties such as porosity and permeability, other important heterogeneities can be difficult to include, such as the possibility of enhanced flow in fractures or along fault zones.Likewise, anecdotal information regarding, for example, the presence of thief zones during waterflooding, are not often explicitly incorporated into the reservoir flow model.Similarly, for compacting reservoirs, the changes in permeability that accompany deformation, and the likely introduction or exacerbation of permeability anisotropy is, if included at all, generally done so on an ad hoc basis, rather than being based on laboratory observations or other more sophisticated modeling methods to determine flow properties.These are limitations that need to be recognized, and in some cases, examined in a quantitative fashion.
Unfortunately, even when the geomechanics problem is well-defined, and even when there is reasonable experimental rock mechanics data, there are many sources of uncertainty that can lead to uncertain results.Uncertainties abound in nearly every aspect of a well-posed geomechanics problem.These uncertainties can be found in physical properties, constitutive laws relating stress to strain, their associated material parameters, initial and boundary conditions, the many types of static or dynamic loads encountered in a geomechanics setting, and even in the construction of the finite element mesh from geologic or other numerical data.Many of the uncertain variables can be categorized as random variables and quantified in terms of probability density functions, which may also exhibit temporal or spatial variation caused by natural geologic heterogeneity and anisotropy.Others can only be given in terms of probable ranges of values.Given these many uncertainties, the geomechanics analyst potentially faces the task of calculating response values and in assigning some level of confidence to these response values.It is sometimes convenient to categorize uncertainties into those that are irreducible and those that are reducible.Irreducible uncertainties are intrinsic to nature and beyond one's ability to control.Irreducible uncertainties are thus caused by inherent variability.In contrast, reducible uncertainties are caused by estimation errors, modeling error, and human error.While these errors cannot be completely avoided, additional testing, more accurate modeling, and better quality control measures can attenuate their effects to a large extent.
There has been a resurgence of interest in the use of probabilistic structural analysis methods in more conventional geotechnical engineering applications (e.g., Kulhawy and Prakoso, 2000) prompted by an increased awareness of the need to go beyond customary deterministic procedures, based on a simple go/no-go basis.Because deterministic geomechanics analysis tools such as JAS3D have reached a highly developed state of maturity and because computer systems have entered an era of truly high-performance computing, this should motivate application of these techniques to geomechanical modeling.The issue of reliability will be an increasing concern in petroleum geomechanics when the cost of failures rises to a critical level.A case in point is the deepwater GoM oil fields where drilling and completion costs for a single field development can run into billions of dollars such that unexpected well failures can profoundly alter the economics of the development.Real-world needs call for the ability to go beyond traditional deterministic methods.Today's well casing designer and reservoir manager must be concerned with the ability to predict risk, minimize risk within limits of cost and need, and understand which features of a particular design lead to risk.Development of procedures for answering those questions efficiently and accurately by probabilistic means should be a goal of the industry when working in challenging environments.

Implementation of Far-Field Stresses and Stress Initialization
Geomechanics problems require special procedures to handle initial in situ stresses.In general, in the normal direction, either displacement or stress can be applied to the mesh, but not both, as this results in an ill-posed problem.Since it is desirable to use zero-displacement boundary conditions to model infinite extent, this requires a special procedure to establish an initial stress field that exists at zero displacement such that field applications can be modeled realistically.The JAS3D code allows self-equilibrating initial stress fields to be specified for each finite element, or uses an equilibrium iteration procedure to achieve equilibrium if the initial stress field is not in equilibrium (which it rarely, if ever, is for real applications).The resulting displacements are then set to zero for subsequent analyses.In general, because of the nonlinearity of the material models, the assumed initial stress state can affect the ultimate numerical calculations of stress and strain.For example, for material models with a state variable, such as the continuous surface plasticity model, the yield surface is subjected to differing hardening depending on the assumed initial stress state.Assessment of the importance of the potential model sensitivity can be established by parametric studies.

Integration of Geomechanical Modeling Results with Other Analysis Tools
The usefulness of geomechanics finite element simulations is enhanced by the ability to port the analysis results directly into other software packages so that one can visualize and examine the geomechanics results in the context of other data of interest for field development and future reservoir management.For example, it can be beneficial to examine potential well paths in three dimensional visualizations of shear stress, or to examine projected shearing displacements in the overburden along well paths above a compacting reservoir, or to predict the stress state better when entering or exiting salt bodies that overlie a reservoir.While we have developed a translator to port analyses from the EXODUSII database format (Schoof and Yarberry, 1995) used in the Sandia finite element codes into General Mesh Viewer (GMV), a public domain highly flexible three-dimensional visualization software package developed at Los Alamos National Laboratory (Ortega, 2001), it is also equally feasible to develop translators for other visualization software used by the industry such as : -IBM Open Visualization Data Explorer (http://www.research.ibm.com/dx/) or -GoCad (http://www.ensg.u-nancy.fr/GOCAD/).

CONCLUSION
The essential question underlying much of the above discussion is: how much is enough?This underpins virtually every area that is discussed.For example, when is it necessary to develop three-dimensional finite element meshes as opposed to simpler, and smaller, two-dimensional meshes?Ideally, at an early stage in the planning of the modeling effort, geologists are included in the discussion so that the numerical modeler can capture the requisite detail.For some structural geomechanical calculations such as at WIPP, where one may be interested in the closure of highly elongated disposal rooms, it may be adequate to consider a plane strain idealization.However, analyses of the behavior at entry-entry intersections naturally require three-dimensional analyses.For geomechanical reservoir modeling, it is more difficult to imagine capturing the necessary structural features in a two-dimensional model.For example, even in cases where the geologic structure is more or less symmetrical, such as at the diatomite fields, if there is a high degree of well interaction, then it again becomes critical to base decisions on three-dimensional modeling.On the other hand, for some near-field applications, such as the casing design analyses presented above, two-dimensional calculations that assume plane strain should be adequate.
These issues are less clear-cut when making determinations regarding the constitutive models to be developed for both the producing as well as overburden formations, or the necessity of including structural features that may be important such as faults, joints, or weak bedding planes, or the need to update the reservoir flow properties due to deformation of the reservoir formation.Unfortunately, the answers here depend quite critically on the questions that must be answered, and on how precisely those answers need to be known.To deal with these issues, analysis of important uncertainties is required.For example, in some cases, the only way to know that it may be reasonable to ignore the presence of identified fault planes is by constructing a preliminary model in which the role and behavior of the fault planes may be determined.Likewise, the influence of the assumed initial stress state on the model results may often best be determined by direct parametric studies.The realworld needs of developing difficult reservoirs, where the cost of failure can be very high, call for the ability to go beyond traditional deterministic methods and enter the realm of probabilistic analysis.
While geomechanical modeling has evolved in its sophistication and application, the above limitations need to be recognized.Fundamental to this is the concept that the numerical model is an approximation of reality.The modeling and simulation tools will always provide an answer; however, there is no guarantee that that answer is reasonable.Therefore, benchmarking and validation of the numerical modeling should be integral parts of the overall effort.Likewise, the importance of experience on the part of the analyst cannot be underestimated.
Figure 1 Schematic representation in I 1 (equal to three times the mean stress) -√J 2 (equal to the second deviator stress invariant) space showing (a) the shear yield surface and hardening cap surfaces of the traditional cap plasticity model, and (b) subsequent hardening surfaces of the (new) continuous surface plasticity model.

Figure 2
Figure 2 Figure 3Measured and predicted volumetric strains as a function of axial stress for a series of triaxial compression tests performed on diatomaceous sand from the Lost Hills oil field, California, for the (new) continuous surface plasticity model.The continuous surface plasticity model predicts the transition from compaction to dilation that occurs prior to peak stress (1 psi = 6894.8Pa).

Figure 4
Figure 4 Schematic representation of surface deformation and shear localization accompanying fluid extraction from a reservoir (after Segall, 1989).

Figure 5 Strain
Figure 5Strain versus time for WIPP and SPR salts at a stress difference of 10 MPa (1450 psi) and temperature of 240°F.The WIPP salt is a bedded salt formation in southern New Mexico, and the SPR salts are from Big Hill, Jennings, Avery Island, Weeks Island, Bryan Mound, West Hackberry, Bayou Choctaw, and Moss Bluff, all salt domes located along the US Gulf Coast.

Figure 11 Finite
Figure 11 Finite element mesh constructed for geomechanical simulations of Sections 4, 5, 32, and 33 at the Lost Hills Oil field, California.The top of the model is at the earth's surface, and the lines above and to the right side indicate the areal and vertical extent of production, respectively.The model was meshed from marker data and the spatial dimensions are 1975 × 15 000 feet × 5400 ft in depth.The mesh contains 258 962 elements and 300 000 nodes.The model also contains contact surfaces that are located along two bedding planes that the field data indicate serve to localize casing damage.(1 ft = 0.3048 m). Figure10 Figure 12Nodal displacements at years 4, 10, 11, 12, and 17 as a function of depth in the horizontal x-direction (left), horizontal y-direction (middle), and vertical z-direction (right) corresponding to a well location for the Section 33, Belridge Field simulation.The horizontal lines bound the reservoir formation, with the upper line also being a frictional contact surface.Two additional contact surfaces (not indicated) are located farther up in the overburden (1 ft = 0.3048 m).
Figure 14 Sequence showing (a) the finite element mesh and expanded views for (b) uncemented elliptical borehole with casing already set and (c) cemented elliptical borehole with casing and cement already set; (d-f) von Mises stress in the uncemented casing at (d) initial contact of salt with casing; (e) initial yielding in casing; (f) complete encapsulation of casing by salt; and factor of safety (defined as strength divided by von Mises stress) in (g) casing and (h) cement at end of service lifetime for the case of the cemented hole.