Iterative Coupled Analysis of Geomechanics and Fluid Flow for Rock Compaction in Reservoir Simulation

Iterative Coupled Analysis of Geomechanics and Fluid Flow for Rock Compaction in Reservoir Simulation — Conventional reservoir simulators calculate the effect of rock compaction on pore volume change through the concept of rock compressibility under a defined loading condition (hydrostatic or uniaxial strain). This approach usually is appropriate for reservoirs with competent rock. For weaker formations and complicated rock compaction behavior, however, a coupled analysis of geomechanics and multiphase fluid flow may be required for obtaining more rigorous and accurate solutions from reservoir simulation. In general, computational efficiency and convergence of numerical solutions are two critical factors in order to make coupled analysis economically and numerically feasible for practical field applications. Oil & Gas Science and Technology – Rev. IFP, Vol. 57 (2002), No. 5


INTRODUCTION
Conventional reservoir simulators calculate the pore volume change through the rock compressibility concept with the assumptions that in the reservoir, the total stress (i.e., the effective stress + fluid pressure) is constant; the loading condition is identical to the laboratory condition under which the rock compressibility was measured; and the local bridging effect around a gridblock from its neighboring gridblocks can be ignored. Also, conventional reservoir simulators ignore the interaction effects between the reservoir and its surrounding regions such as overburden, underburden, and sideburden. The rock compressibility approach is generally applicable to reservoirs with competent rock, laterally uniform rock properties, and reservoir properties such as permeability, which are insensitive to the change of stress state. These assumptions, however, may restrict conventional reservoir simulators from analyzing reservoirs of complex geomechanical behavior. A coupled analysis of reservoir simulation and geomechanics has significant advantages over conventional reservoir simulation for these cases. With the integration of geomechanics, reservoir rock of complicated constitutive behavior can be rigorously simulated without any geomechanical limitation. The geomechanics part of the coupled analysis allows one to accurately calculate reservoir rock deformation and pore volume change. Through geomechanics, the coupled analysis is capable of handling the effects of stress-sensitive properties and heterogeneity in the reservoir. The coupled analysis also takes into account the effects of surrounding regions/reservoir interaction and local bridging in the reservoir. Therefore, for geomechanical conditions that conventional reservoir simulation cannot handle, coupled analysis of reservoir simulation and geomechanics should be applied.
There are many cases that conventional reservoir simulators may not be able to handle by using the rock compressibility concept and reservoir property change that is related to the stress change. Typical cases are reservoirs with weak rock such as chalk and unconsolidated sand. Weak rock usually exhibits complicated constitutive behavior and results in high reservoir compaction and stress arching in the overburden. It is necessary to employ a coupled analysis for accurately quantifying pore volume change and reservoir compaction. Another case is for reservoirs with high degrees of heterogeneity in rock mechanical properties, where the effect of local bridging may become pronounced in the reservoir. It requires the geomechanics part of a coupled analysis to rigorously compute pore volume change under a non-uniform compaction condition. Similarly, for reservoirs with stress-sensitive properties such as permeability, permeability change can be related to the stress change that is calculated by the geomechanics model.
In general, the coupled analysis of geomechanics and reservoir fluid flow can be classified into three types (e.g., [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]). They are loosely coupled, iterative fully coupled, and simultaneously fully coupled. The simultaneously fully coupled analysis solves the coupled system of equations of geomechanics and reservoir multiphase flow simultaneously. The resulting solution is self-consistent and considered to be the true solution of the coupled problem. However, the primary disadvantage of this approach is that its computational cost is extremely high compared to other types of coupled analyses. Also, it may be difficult to develop a robust and effective numerical procedure for simultaneously solving the coupled system of equations that combines geomechanics and reservoir multiphase flow. On the other hand, depending upon the degree of coupling that varies between one-way coupling to explicit coupling (e.g., [1][2][3][4]), the loosely coupled analysis has low computing cost and at best provides only an approximate solution to the true solution of the coupled problem. The iteratively coupled analysis decomposes the coupled system of equations into two subsystems of equations that correspond to the governing equations of the reservoir simulator and of the geomechanics model. It solves the coupled system iteratively by exchanging shared state variable values between the reservoir simulator and the geomechanics model through a data exchange interface, while the reservoir model and the geomechanics model separately solve their corresponding subsystems of equations. When the iterative process converges, the iterative fully coupled analysis also reaches the true solution of the coupled problem that integrates geomechanics and reservoir multiphase flow.
Though both iteratively coupled and simultaneously coupled analyses obtain the true solution of the coupled problem, the iterative fully coupled analysis has several advantages over the simultaneously fully coupled analysis. First of all, the iteratively coupled analysis has better 486 In this paper, an iterative procedure for coupled analysis of geomechanics and multi-phase flow in reservoir simulation is proposed for large-scale, full-field, 3D problems. The proposed procedure is general and effective for handling reservoir rock with complicated constitutive behavior of rock compaction and permeability change as well as for simulating various reservoir production scenarios. Descriptions of model formulations, constitutive equations, solution procedures, and strategies for enhancement of computational efficiency are presented in the paper. To demonstrate the capability of the developed procedure for iterative coupled analysis, several problems including a large-scale field example were studied and are presented.
computational efficiency than the simultaneously coupled analysis. In addition, through the data exchange interface, the iteratively coupled procedure can easily be implemented with existing reservoir simulators and existing geomechanics models. A number of sophisticated reservoir simulators and geomechanics models are readily available for implementation of the iterative procedure. On the contrary, it will generally take great efforts to develop and implement a simultaneously coupled procedure. Finally, new computing technologies and numerical methods can be independently developed and implemented for the reservoir simulator and the geomechanics model used in the iterative coupled analysis. For example, the reservoir simulator and the geomechanics model typically use the finite difference method and the finite element method for solving their governing equations, respectively. They also use different linear solvers that are best to solve their corresponding systems of linear equations. Parallel computing technology can be independently developed for the geomechanics model and the reservoir simulator. Consequently, this advantage generally results in a very robust and effective computing procedure for the iteratively coupled analysis. Based on the advantages described above, we have developed a computational procedure for iterative coupling in this paper.
Until now, most of the fully coupled models proposed were used for analyzing small-to medium-scale field problems (e.g., [7-9, 15, 16, 19]). Employing the fully coupled model of reservoir simulation and geomechanics for full-field, 3D field problems was generally considered impractical. Because of the difference in governing equations and numerical methods, geomechanical modeling is normally much more computationally expensive than reservoir simulation. In addition, for the coupled problem, the mesh size of the geomechanical model is larger than that of the reservoir model, since the reservoir model grid is only a part of the geomechanical model grid (i.e., the reservoir section), while the geomechanical model grid includes additional regions such as overburden, underburden, and sideburden. As a result, the computational cost for solving large-scale field problems using the fully coupled model of reservoir fluid flow and geomechanics tends to become very high. Advances in numerical simulation techniques and high-performance, lowcost computing technologies, however, make it economically feasible to conduct fully coupled analyses that integrate geomechanics and reservoir simulation. By utilizing these advanced techniques and computing technologies, the iterative, fully coupled procedure proposed here can efficiently analyze field-scale, 3D problems that integrate reservoir simulation with geomechanics.

ITERATIVELY COUPLED ANALYSIS
In this paper, we propose an iterative, fully coupled procedure that integrates reservoir simulation with geomechanics in a generalized fashion. Two computer models, a reservoir multiphase flow model and a geomechanics model, are linked together through discretized coupled equations. These discretized coupled equations contain state variables shared by both reservoir and geomechanics models. These shared variables are porosity and fluid pressure for each gridblock. Within each time increment, both models iteratively provide each other numerical values of needed state variables on a gridblockby-gridblock (element-by-element) basis. The reservoir section in the geomechanics model has the same mesh configuration as the reservoir model. The state variable value needed by the geomechanics model and provided by the reservoir model is pressure. The state variable value needed by the reservoir model and given by the geomechanics model is the porosity defined by constant bulk-volume, which is traditionally used in reservoir simulators with fixed gridblock volumes and in this study, is called the reservoir-simulation porosity to distinguish it from the true porosity. The reservoir model also needs the derivative of porosity with respect to pressure, which can be calculated from the geomechanics model. This iterative process of performing reservoir simulation and geomechanics computation by the two separate models and passing needed state variable values between the two models continues until convergence criteria are satisfied. Following convergence, the coupled analysis marches to the next time step and repeats the iterative process described above.
A flowchart illustrating the iterative, coupled procedure is shown in Figure 1. The procedure consists of an initialization phase and a solution phase. At time = 0, both reservoir and geomechanics models are set at their initial conditions. Based on the initial fluid pressure distribution provided by the reservoir model, the geomechanics model solves the equation of equilibrium to establish the initial stress state of the modeled system that is at equilibrium under the prescribed initial conditions. Then, the geomechanics model gives initial porosity values and initial derivatives of porosity with respect to pressure to the reservoir model. At this stage, the initialization phase is complete and the reservoir model is ready to start the first time-step calculation.
In the solution phase, within any given time increment (i.e., dt = t n -t n-1 , n = 1, 2, 3, etc.) the reservoir simulator and the geomechanics model solve their nonlinear systems of equations separately and provide needed state variable values to each other through a coupled iteration loop with an iteration number counter k (k = 1, 2, 3, etc.). Both reservoir and geomechanics models solve their nonlinear systems of equations by the Newton method. However, in this proposed procedure, we make the Newton iteration loop of the reservoir model coincide with the coupled iteration loop by sharing the same iteration counter. As a result, for coupled iteration k, the reservoir simulator solves its nonlinear system at iteration k for its Newton loop, while the geomechanics model solves its nonlinear system until its Newton loop reaches convergence. The coupled iterative process continues until the Newton loop of the reservoir simulator meets its convergence criteria and converges. As a consequence of sharing the same iteration counter k with the Newton loop of the reservoir simulator, the coupled iteration loop is also converged because the geomechanics model always reaches its converged solution for every iteration k. After the solution update, the iteration counter is reset to 1. If time is less than the maximum time (t max ), the coupled analysis marches to the next time increment and repeats the same iterative process. Otherwise, the solution phase for the iteratively coupled analysis is stopped.
In the following, descriptions of the main components (the reservoir model, the geomechanics model, and the constitutive model of reservoir rock) in the iterative, coupled procedure are presented.

Governing Equations
The governing equations for the classical black-oil reservoir model are presented here for simplicity (e.g., see [20]), however, the coupling analysis described in the paper is general and can be applied to any reservoir model formulation. Based on the conservation of mass for each phase and Darcy's Law, the following multiphase flow equations are obtained. (1) where: λ l is the transmissibility of phase l and w, o, and g refer to the water, oil, and gas phases, respectively; B l = formation volume factor of phase l; p l = pressure of phase l; q l = volume of stock tank phase l produced per unit of reservoir volume per unit time; R s = solution gas-oil ratio; S l = saturation of phase l; t = time; z = depth measured positive downward; γ l = weight density of phase l; φ = reservoir-simulation porosity.
Additional relations required to complete the formulation of governing equations are: where f 1 and f 2 are functions for the water-oil and gas-oil capillary pressures, respectively.
The governing equations for the geomechanics model are derived based on the balance of linear momentum and the effective stress law [21]. The equations are stated as follows: Initialization phase t 0 = 0 Reservoir simulator Geomechanics model Reservoir simulator Geomechanics model Newton loop converged?
Flowchart of the iteratively coupled analysis.
In order to derive the full set of governing equations, Equation (7) must be supplemented with an appropriate constitutive equation for the solid porous matrix. In general, the constitutive equation is in the form: = solid strain; ε . = solid strain rate; Κ denotes a set of history dependent tensorial quantities. Assuming the solid grain is incompressible, the solid strain rate is related to solid velocities as: (9) where ∇v s 0 = symmetric part of the solid velocity gradient. These nonlinear constitutive equations, Equations (8) and (9), are usually appropriate to describe the nonlinear behavior of weak reservoir rocks. With the assumption of small strain in a time increment, the solid strain tensor is related to solid displacement as: (10) where ∇u 0 is the symmetric part of the solid displacement gradient and u is the solid displacement vector. Introducing Equations (8, 9, and 10) into Equation (7), we get a partial differential equation with the solid displacement vector u as the unknown variable. After solving the partial differential equation for u, σ' and ε can be calculated from Equations (8) and (10), respectively.
Based upon the compaction behavior measured in the laboratory and observed in the field, the constitutive model required to properly model reservoir rock compaction could vary from a simple linear-elastic model to a sophisticated elasto-plastic model. Typical well-structured 3D, nonlinear geomechanics models can easily implement various constitutive models for simulating different types of rock compaction behavior. In this study, we used a constitutive model that is capable of simulating the elasto-plastic compaction behavior to exemplify the generality of the fully coupled procedure presented here. The nonlinear crossanisotropic constitutive model developed is based on a hypoelastic/hypoplastic, 3D formulation [22] that is capable of reproducing uni-axial strain test data over a wide range of loading conditions. Thus, typical nonlinear elasto-plastic behavior of weak rocks can be phenomenologically simulated by the proposed model. Application examples of using the proposed constitutive model for rock compaction will be presented later in the paper.

Implementation of the Iterative Coupled Procedure
The iterative, coupled procedure solves the governing equations of the reservoir model and the geomechanics model separately and iteratively, as illustrated by the flowchart in Figure 1. Numerical solutions for Equations (1)(2)(3)(4)(5)(6) and for Equations (7-10) are straightforward because numerical simulation techniques for both petroleum reservoir simulation and computational geomechanics have been highly developed. The primary unknowns, solved from Equations (1-6) together with the values of φ and ∂φ/∂p provided from the geomechanics model, are p o , p w , p g , S o , S w , and S g . The primary results, solved from Equations (7-10) together with the values of p provided from the reservoir simulator, are u, σ', and ε.
Finite difference methods such as implicit + pressureexplicit saturation (IMPES) and fully implicit simultaneous solution methods [20] are used for solving the reservoir flow Equations (1-6) to obtain pressure and saturation distributions in the reservoir. The finite element method along standard lines [23] is used to solve Equations (7-10) for the displacement, stress, and strain fields in the computational domain of the geomechanics model, which includes the reservoir and additional regions (the overburden, sideburden and underburden).
Within a given time increment, for each coupled iteration through a data exchange interface, pressure values obtained from the reservoir model are directly passed to the geomechanics model and porosity values and derivatives are passed from the geomechanics model to the reservoir model. The data exchange interface uses uniquely named disk files. After each data file is written and closed, a second uniquely named file containing time, time step number, and iteration number is written. The existence of the second, small file signals that the model awaiting data may open and read the new data file. It should be noted that the mathematical system of the reservoir model is based on the Eulerian description with a fixed mesh configuration that is timeindependent, while the mathematical system of the geomechanics model is based on the Lagrangian description with a deformable mesh configuration that is changing with time. Thus, the true porosity calculated from the geomechanics model cannot be directly passed to the reservoir model. In order to maintain numerical consistency in the definition of porosity used in the reservoir model, the true porosity of a given element has to be converted to the reservoir-simulation porosity defined by the constant bulk volume for the corresponding grid block in the reservoir model. The equations that relate the true porosity and the reservoir simulation porosity to the volumetric strain of a given element in the geomechanics model are: where φ 0 is the initial porosity at time = 0; ε v = tr ε = volumetric strain of an element, the prefix tr denotes trace. At time = 0, the initial reservoir-simulation porosity and the initial true porosity are the same. Note that ε v is negative for compaction. Using Equation (12) the reservoirsimulation porosity for each grid block was obtained from the geomechanics model and passed to the reservoir model. Substituting Equations (8) to (10) into Equation (7), the finite element discretization of Equation (7) was expressed in terms of the nodal displacement vector field u using the Galerkin method. The unknowns of the field equations are related to their nodal values by: where N is the shape function, and the finite element procedure yields the following discretized form of the field equations: where n s are nonlinear vector-valued functions of the nodal solid displacement vectors; f p is a vector arising from the pressure gradient term in Equation (7) as: is the elemental contribution to node a (a = 1, 2,.., nn) from direction i (i = 1, 2, and 3) and element e through the weak form formulation; nn = total number of nodes; Ω e = domain of the finite element e; and f s is the term contributed from boundary conditions through the weak-form formulation.
Since both reservoir and geomechanics models use the same grid in the reservoir, fluid pressure values obtained from the reservoir model for a given time increment, were used directly in Equation (15) for all gridblocks (i.e., elements). Note that p in Equation (15) is at the cell value from the reservoir model. The solution of the displacement from the finite element procedure can be used to compute ε from Equation (10). The volumetric strain ε v can, in turn, be calculated from the strain tensor ε. Then, the geomechanics model uses the volumetric strain to calculate the reservoir porosity through Equation (12) for all gridblocks and provides the porosity values to the reservoir model.

APPLICATION EXAMPLES
Two examples are presented to illustrate the performance of the coupled reservoir and geomechanics models. The first case presents the depletion of a five-spot area of a field followed by waterflooding. The second example is a large full-field model that includes production, gas injection, and water injection. In both of these examples permeability has been assumed to not be a function of stress. For cases where stress-dependent permeability is important, this feature can be included by entering tables of permeability versus porosity or by using functions that relate permeability to stress (e.g., see [19]). Comparisons between coupled and uncoupled results for specific examples were given in Reference [19].

Five Spot Example
This example illustrates the coupled simulation of a five-spot area of a layered reservoir. The grid for the reservoir model was 11 × 11 × 5 and the grid for the geomechanical model was 11 × 11 × 15 with eight layers in the overburden and two in the underburden. Overburden and underburden layer thickness was 381 m. The areal grid (11 × 11) in the x-and ydirection was equally divided and the gridblock dimension in both directions was 61 m (i.e., an area of 671 m × 671 m).
Depth to the top of the reservoir was 3048 m and depth to the bottom of the underburden was 3880 m. The field was produced under primary depletion for the first four years of operation and then waterflooding was initiated. The maximum rate of the producing well, which is located at the center of the five-spot was set equal to 397.5 stock tank m 3 /day and the minimum flowing bottomhole pressure was set equal to 6.895 MPa. The maximum injection rate for each of the four water injectors  associated with this five-spot was set at 238.5 m 3 /day with a maximum bottomhole injection constraint of 51.71 MPa. Oil and water compressibility were 1.45 GPa -1 and 0.508 GPa -1 , respectively. Oil and water viscosity were 0.5 and 0.35 cP, respectively. Vertical to horizontal permeability ratio was 0.1. The connate water and residual oil saturations were 0.1 and 0.3, respectively. Normalized saturation squared curves were used for relative permeability curves. Other pertinent reservoir data for this example are presented in Table 1.
Oil rate and water cut versus time are presented in Figures 2a and 2b, respectively. Oil production rate is constant for the first 10.5 years before going on decline at approximately the same time as water breakthrough, as shown on Figure 2b. Water injection rate for the five-spot versus time is shown on Figure 3a. Figure 3b shows that calculated reservoir pressure in the top well cell declines from an initial pressure of 48.26 MPa to 16.34 MPa during depletion and then increases to 34.47 MPa during the waterflood. Figure 4 presents calculated total reservoir compaction at the central location as a function time. The plot clearly indicates that total reservoir compaction increases with the decrease in hydrocarbon pressure under depletion for the first four years of operation. During the period of waterflooding, the elastic recovery from the increase in pressure was negligible and thus, permanent deformation took place in the reservoir. This result verifies that the constitutive model of reservoir rock, used in the geomechanics model, can properly simulate the hysteresis effect caused by plastic deformation. Oil rate as a function of time for the five-spot example.

Figure 2b
Water cut as a function of time for the five-spot example. Hydrocarbon pressure (MPa)

Figure 3a
Water injection rate as a function of time for the five-spot example.

Figure 3b
Calculated hydrocarbon pressure as a function of time for the five-spot example.

Field Example
A large anti-cline field with a history of compaction and subsidence was selected for this example. Figure 5 shows oil production history from the field including 41 active production wells. A three-component, compositional model [24]  Drive mechanisms in the early operation of the field were oil expansion, limited water influx, compaction drive, and solution gas drive as the reservoir pressure fell below the bubble point. Gas injection in eight wells was initiated after 4 years of production and continued to some extent to the end of the 17-year run, as given in Figure 6a. Limited water injection in four wells was also included in parts of the model during the last 7 years of the simulation, as presented in Figure 6b. Figure 7 shows calculated average reservoir 492 Figure 4 Calculated total reservoir compaction versus time for the five-spot example.

Figure 5
Oil rate versus time of the field example.

Figure 6a
Gas injection rate versus time of the field example.

Figure 6b
Water injection rate versus time of the field example.  Figures 8 and 9, respectively. Some simulation results from the geomechanics model are given in Figures 10 through 12. Figure 10 shows subsidence at the central location of the field versus time. Figures 11a  and 11b show contour maps of seafloor subsidence and total reservoir compaction at year 17, respectively. Figure 11b indicates that total compaction contours are not uniformly distributed around the reservoir center and there exists some local maximum compaction areas, because of non-uniform distributions in initial porosity, mechanical properties of reservoir rock, and pore pressure. The overburden tends to smooth the downward movement induced by reservoir compaction and gives smooth seafloor subsidence contours around the reservoir center, as illustrated in Figure 11a. Vertical stress change in the overburden layer eight at year 17 is shown in Figure 12a. The contours shown in Figure 12a verifies the stress arching in the overburden with a decrease in vertical stress in the central area and an increase in vertical stress around the flank region. Effective vertical stress change in the reservoir layer 5 at year 17 is presented in Figure 12b. Again, the contours on Figure 12b indicate that effective vertical stress change in the reservoir is not uniformly distributed and there exists local areas of maximum and minimum stress changes.   Gas-oil ratio calculated by the couple analysis for the field example. Water cut calculated by the coupled analysis for the field example.      This large-scale, 3D field example has a modeled regime that covers an area of 19.17 km × 18.38 km with a depth of 4.57 km. The reservoir model has 38 458 (67 × 41 × 14) gridblocks. The geomechanics model has about 200 000 degrees of freedom and 65 928 (67 × 41 × 24) elements that consist of the reservoir section and the overburden, sideburden and underburden regions. Note that the 3D mesh was not uniformly divided in the areal grid as well as in the depth dimension for both models. The mesh configuration and geological descriptions of the reservoir section in the geomechanics model (i.e., 38 458 elements) are identical to those of the reservoir model. The reservoir simulated in the study has high degrees of heterogeneity in initial porosity distribution and mechanical properties. To verify the robustness of the proposed procedure, reservoir rock with complicated compaction behavior was studied in the geomechanics model. The stress-strain relationship of the reservoir rock is a function of initial porosity and was modeled by a hypoelastic/hypoplastic constitutive model [22] for each porosity value. This field example took 181 time steps and 235 coupled iterations (i.e., 235 Newton iterations in the reservoir model) for the 6024-day simulation. Thus, the iterative coupled procedure performed the simulation run for the field example in a very robust and effective manner.
The total elapsed time for the entire run took only a few hours. These results indicate that the proposed computational procedure can be practically used for analyzing fieldscale problems that integrate reservoir simulation with geomechanics. Figure 13 Schematic of the computation process for the iteratively coupled analysis.

DISCUSSION
Simulation runs for application examples presented in this paper were conducted on a cluster of Pentium personal computers (i.e., nodes). As illustrated in Figure 13a, two nodes on the cluster were used for the coupled analysis with the reservoir simulator running on node 1 and the geomechanics model running on node 2. The two-headed arrow indicates the data exchange and intercommunication between the two nodes each coupled iteration. Though the reservoir simulator and the geomechanics model are running independently on their own nodes, the reservoir simulation and geomechanical computation are sequentially performed (see Fig. 1) each coupled iteration and are coordinated by a hand-shacking protocol. In general, the computing speed of the geomechanics model is considerably less than the reservoir flow model. Based on the results of application examples, the ratio of reservoir to geomechanics model run times was less than 0.05. As described in the Introduction, reasons for this disparity include the increased number of gridblocks required in the geomechanics model to simulate the overburden and underburden in addition to the reservoir; the finite element discretization of the primary equations of the geomechanics model; and the stress integration of the constitutive equations for plasticity. Also, the geomechanics model simultaneously solves three times more equations per gridblock than the reservoir flow model when it is run in IMPES mode. The geomechanics model performs the simultaneous solution for deformation in all three directions ( a 3-component vector), while the reservoir flow model simultaneously solves the system of equations for pressure that is a scalar. As a consequence, the geomechanics model is the bottleneck in computation for the iterative, fully coupled analysis. Any improvement in computing speed of the geomechanics model can result in increased computational efficiency for the iteratively coupled analysis.
There are two approaches that can drastically increase the computing speed of a geomechanics model. One is to implement the fastest linear solver available for the geomechanics model since 70 to 80 percent of the run time is spent in solving the system of linear equations. We have worked on this approach in the past 6 years. Depending upon the types of geomechanical problems being solved, a speedup by a factor of 2 to 5 was achieved. This approach, however, will eventually reach a limit as the development of faster linear solvers becomes more challenging and difficult. The other approach is to implement parallel computing technology for the geomechanics model to perform geomechanical calculations in parallel. Parallel computing on a cluster of personal computers (PCs) interconnected by a private high-speed network is an economical and proven technology for high-performance computing. The performance of parallel computing on a cluster generally increases with the number of PCs (i.e., nodes) in the network. Parallel geomechanics model Therefore, a multi-fold increase in computing speed is easily achievable. The limitation of parallel computing performance is normally constrained by the communication among the nodes. Work on development and implementation of parallel computing for the geomechanics model was started. Along the standard lines of parallel computing on a cluster of PCs, the domain decomposition method with the message-passing interface (MPI) [25] was employed for developing the parallel geomechanics model. The parallel geomechanics model was implemented with a parallel iterative linear solver that is based on the preconditioned conjugate gradient algorithm. The parallel geomechanics model was verified by a number of test problems. In general, the expected performance in speedup was obtained. The preliminary results for solving large-scale problems appear quite promising. For example, a field-scale, 3D, geomechanical problem was studied by the parallel computing method developed. The mesh configuration of the problem is a brickshaped regime that is composed of 50 × 50 × 25 (62 500) hexahedral elements. The problem has a single point-load, which is a function of time, applied on the top surface of a corner element and the material modeled in the problem is treated as linear elastic. The performance of parallel computing for this problem is shown in Table 2. Speedup is the ratio of the run times between the parallel run and the single-node run. Efficiency is defined as the speedup in terms of the percentage of the number of nodes used. The results clearly show that computing speed increase by an order of magnitude was accomplished when 16 nodes on the cluster were used. Work is in progress to integrate the iterative, fully coupled procedure with the parallel geomechanics model. It is expected to further significantly increase the computational efficiency of the iteratively coupled analysis. The schematic of the computation process on a cluster of PCs for the iteratively coupled analysis integrated with a parallel geomechanics model is illustrated in Figure 13b. The sequential reservoir simulator is running on a single node, while the parallel geomechanics model is running on multiple nodes (e.g., 4 nodes in the circle). Again the twoheaded thin arrow indicates the data exchange and intercommunication between the reservoir simulator and the geomechanics model each coupled iteration. The two-headed block arrows represent the intra-communication among nodes in the parallel computing procedure of the geomechanics model. As shown previously in Table 2, the parallel computing speed of the geomechanics model increases with increasing number of nodes used in the circle. Thus, the computational efficiency of the iteratively coupled analysis will be increased drastically by the proposed scheme (Figure 13b) of sequential reservoir simulation integrated with parallel geomechanics computation.

CONCLUSIONS
• An iterative, fully coupled procedure that integrates reservoir multiphase flow simulation with geomechanics was developed and implemented into a reservoir simulator and a geomechanics model. The iterative fully coupled procedure can rigorously and accurately calculate pore volume change for reservoirs with complex geomechanical behavior and high degrees of heterogeneity. The procedure properly takes into account the interaction effects between the reservoir and its surrounding regions such as overburden, sideburden, and underburden. • The coupled procedure was verified by practical application examples. The simulation results indicate that the procedure can be used to effectively and efficiently analyze field-scale, 3D problems of fully coupled geomechanics and reservoir multiphase flow. • The iteratively coupled procedure is general and efficient and can be implemented with any existing reservoir simulator and any existing geomechanics model through a data exchange interface. • A parallel computing method was developed and implemented for geomechanics models. The study indicates that this method can increase the speed of geomechanical computation by an order of magnitude. It is expected to significantly increase computing efficiency of coupled analysis when the iteratively fully-coupled procedure is integrated with parallel geomechanics models.