Numerical Modeling of Thermal EOR: Comprehensive Coupling of an AMR-Based Model of Thermal Fluid Flow and Geomechanics

Résumé — Modélisation numérique d’EOR thermique : couplage complet entre un modèle d’écoulement thermique basé sur une discrétisation adaptative et la géomécanique — La modélisation du procédé SAGD ( Steam Assisted Gravity Drainage ) peut impliquer un temps de calcul important lorsque l’écoulement thermique et la géomécanique sont couplés pour tenir compte des variations de perméabilité et de porosité à l’intérieur du réservoir induites par l’évolution des contraintes. Une procédure numérique qui effectue des simulations thermo-hydro-mécaniques, d’une manière efficace, est présentée. Cette procédure repose sur un processus de couplage itératif entre un simulateur réservoir thermique basé sur la méthode des volumes finis et un simulateur géomécanique basé sur une discrétisation par éléments finis. Une caractéristique forte de cette procédure est qu’elle permet de traiter des cas où les simulations de réservoir sont réalisées en utilisant un raffinement de maillage adaptatif. Elle fournit ainsi une description précise de l’évolution du front de vapeur d’eau et permet de prendre en compte les effets géomécaniques sans effectuer les simulations géomécaniques sur un maillage raffiné. L’efficacité de cette procédure de couplage est illustrée sur un cas SAGD synthétique mais réaliste. Abstract — Numerical Modeling of Thermal EOR: Comprehensive Coupling of an AMR-Based Model of Thermal Fluid Flow and Geomechanics —


INTRODUCTION
The production of heavy oils and bitumen by injection of steam involves temperature and pressure changes within the reservoir which may modify the stress state of the porous medium. As a result, the porosity and permeability fields can change during the production period. Several works point out this phenomenon. For instance, the increase of permeability in oil sands through a shear-stress modification has been observed by means of laboratory experiments (Touhidi-Baghini, 1998). On the other hand, numerical simulations of SAGD operations underline the influence of geomechanical phenomena on the oil recovery (Ito and Suzuki, 1996;Collins et al., 2002;Lerat et al., 2009).
These observations induce more and more reservoir engineers to couple geomechanical and fluid-flow models when dealing with thermal processes. Different approaches have been proposed for this coupling. In the fully coupled method, the equations of both problems are solved simultaneously within the same simulator in order to obtain consistent solutions. But the integration of geomechanics in fluid-flow models often results in simplifications in either of these two domains. Another solution consists in coupling classical mechanical and reservoir simulators in an external way to be able to use more advanced functionalities, available in both software (Settari and Mourits, 1998;Longuemare et al., 2002;Jeannin et al., 2005;Tran et al., 2008). This second approach has been adopted in the present work.
However, in both cases, taking geomechanical effects into account in reservoir simulations entail long computational times, in particular when the same grid is used for both simulations. For the reservoir simulation, the use of a fine grid is essential in regions where oil flows. Indeed, if the size of the grid blocks is not small enough, temperatures are too averaged and the computed viscosities are too far from their real values to obtain precise production forecasts.
Concerning geomechanics, accurate results can be obtained, in practice, with coarser grid resolutions. This was shown in Guy et al. (2011) on a SAGD test case. In that case, a fine Geomechanical Grid (GG) with 6 175 elements gave the same oil production as a grid composed of 1 235 elements. In addition to a reduction of the computational times of the geomechanical simulation, it was also observed that the optimal time step, used for the coupling of both simulators, is smaller when using a coarser grid. Consequently, the petrophysical properties can be more often updated in the reservoir model without increasing the total simulation time.
On another hand, the use of Adaptive Mesh Refinement (AMR) can also help to reduce the computational times related to the reservoir simulation: indeed, refining dynamically the mesh, only in regions where oil is mobile, strongly increases the efficiency of the simulation (Lacroix et al., 2003;Christensen et al., 2004;Wang et al., 2006). As shown in Mamaghani et al. (2011) on SAGD cases, the reduction of the grid size varies during the simulation and also depends on the criteria used to follow the flow interface. At the beginning of the simulation, when the extent of the steam chamber is small, AMR methods enable to reduce considerably the grid size. The reduction of the cell number decreases afterwards as more oil is heated and flows towards the producer wells but the method still provides substantial savings of computational times. Concerning the following of the flow interface, criteria based, for instance, on saturations give better results than criteria based on temperatures, in particular when high contrasts of permeability exist within the reservoir.
In this work, the approaches proposed by Guy et al. (2011), on one hand and by Lacroix et al. (2003) and Mamaghani et al. (2011), on the other hand, are combined to propose a new coupling scheme of geomechanical and reservoir simulators to obtain precise production forecasts with reduced computational times. The methodology is described in Section 1 and numerical results are presented in Section 2. Conclusions and perspectives are discussed at the end.

METHODOLOGY
The coupling scheme between the geomechanical (Abaqus ™ ) and reservoir (PumaFlow ™ ) simulators is represented in Figure 1.
The simulation starts with an initialization step where the initial states of both the reservoir and geomechanical models are computed to be consistent with respect to the considered equilibriums. The initialization step is performed by considering a Geomechanical Grid (GG) and a Fine reservoir Grid (FG) that has been used to define the models. Then, a first AMR based reservoir grid (AMRG) where the mesh is refined around the well pairs is built from FG as explained later. The reservoir simulator uses the initial porosity and permeability values defined on FG and their upscaled values in the regions where the grid is coarsened.
The production history is divided into periods during which the reservoir and geomechanical simulators run sequentially, iterating several times during one period until the results of both simulations are consistent.
A period i starts with an update of the reservoir grid (AMRG), where the mesh is refined around the well pairs and the flow interfaces and coarsened elsewhere. All data defined on FG like the petrophysical properties (porous volumes, porous volume corrections, permeabilities, endpoints of the relative-permeability curves, etc.) or the current solutions (pressures, saturations, temperatures) to the flow problem are upscaled on the current AMRG to obtain a relevant modeling of the state of the underground at a time t i . A reservoir simulation is then run for period i from time t i to t i+1 . After the simulation, the pressures, saturations and temperatures are transferred from AMRG onto FG and the Geomechanical Grid (GG), using a procedure called the Diffusive Approximation Method (DAM). The geomechanical simulator is then executed to solve the mechanical equilibrium with respect to thermo-poro-mechanics (Biot, 1941;Coussy 1995) and the results are also transferred onto FG using the DAM. Finally, the algorithm tests the convergence of both simulations by checking that: (1) where ϕ g (t i+1 ) stands for the Lagrangian porosity deduced from the geomechanical simulation on FG, ϕ r (t i+1 ) is the Lagrangian porosity used by the reservoir simulator at the end of its run and also transferred on FG, ϕ 0 is the initial Lagrangian porosity and CRIT is the convergence criterion.
If the previous inequality is not satisfied, the permeabilities are updated and the required corrections of porous volume are computed on FG as explained later and a new iterate for period i is made. This loop is performed until convergence is reached. When the convergence is reached, the porosities of the reservoir and geomechanical simulators are consistent and the permeability of the reservoir simulator is also consistent with the geomechanical state of the rocks. Then, the next period is simulated. This procedure is continued until the end of the simulation.
In the proposed methodology, the fields are transferred from a grid to the other one using the diffuse approximation method (Savignat, 2000). For all the reservoir grids, the data are defined at grid centres whereas, for the Geomechanical Grid (GG), they are defined at grid nodes. The diffuse approximation method is then used to find estimates of a scalar field (temperature, pressure, etc.) from a set of nodal values to another one (Nayroles et al., 1991). The starting point is to estimate the Taylor expansion of the studied scalar field at a chosen point by a weighted least-square method which uses only the values at the nearest points. The main advantages of this method are that it only requires sets of discretization nodes and is a local method. It has to be noticed that the diffuse approximation method can be used with various weighting strategies that lead to different and interesting properties. It is also worth noting that the weighting procedure used in the present scheme allows the DAM to interpolate the initial data.
More details about the Adaptive Mesh Refinement (AMR) method used for the reservoir simulation, the geomechanical model are now given in Sections 1.1 and 1.2.

AMR Method for the Reservoir Simulation
The creation of a new AMRG is carried out by a mesh generator at the beginning of a new period. This module also performs the upscaling of the static (porosities, permeabilities, etc.) and dynamic (pressures, temperatures, saturations, etc.) data from FG to AMRG and the creation of the new input files for the simulator.
Note that, once the grid has been updated, it remains fixed during the whole next period. Our AMR method thus differs from the classical approach (see for example Berger and Colella, 1989), which consists in modifying the grid at each time step, according to error criteria. Our approach may be less efficient since the refined zone should be sufficiently large to avoid that the front moves outside this zone during the simulation period. But, on the other hand, it can be applied with any reservoir simulator having a Local Grid Refinement option. It can be seen as a preliminary approach before implementing an AMR option within a simulator.
An AMRG is composed of several levels of refinement denoted by l. These levels range from the finest one, l = 0, which has the same resolution as FG, to the coarsest one l = l max . At a level (l ≥ 1) a grid block is the union of cells located at level l -1. The grouping ratios may differ in each direction. In the presented examples, these ratios are all equal to 3.  Coupling between the geomechanical and reservoir simulators.
An AMRG is made of nested levels where a cell of level l can only be next to a cell of the same level or of level l -1 or l + 1. This property, illustrated in Figure 2, is used to reduce the error and to avoid creating discontinuities in the solutions.
Refined areas are always created in the neighbourhoods of well pairs so that we do not need to upscale the production indexes and we ensure a good precision of the results near the wells. Outside these areas, a fine grid is also used around the flow interfaces.
At the beginning of the simulations (in the pre-heating phase) no-steam is injected, so the only area where the AMRG is refined is around the wells.
The flow interfaces are located using the criteria proposed in Mamaghani et al. (2011). These criteria have been derived from an a posteriori error estimate which has been obtained by assuming that the oil saturation is solution to a simplified transport problem. To our knowledge, rigorous error estimates have not been established until now, even for a thermal dead-oil model. In our applications, heuristics such as saturation or temperature variations are still often used.
The drawback of such heuristics is that their adjustment strongly depends on the resolution of the different grid levels. In Mamaghani et al. (2011), the local estimators, proposed to control the global error of the simplified problem, are proved to be independent on the time step and the mesh resolution. Numerical tests, performed also on SAGD cases, showed that these criteria were indeed less dependent on the grid size.
When the current positions of the flow interfaces are determined, an estimation of their displacement during the period is made from their locations at previous times. The grid is then refined in the regions covered by those displacements. The cells which are not included in the refined zones are coarsened at levels l ≥ 1 according to the nesting rule.
As mentioned before, the update of the AMRG is followed by an upscaling of the static and dynamic data. Different analytic formula can be used for permeabilities. Other data are upscaled with arithmetic means.

Geomechanics
Our approach aims at coupling a reservoir and a geomechanical simulator to update both reservoir and geomechanical permeabilities and set reservoir porosities consistent with geomechanical porosities. The permeability is updated considering the evolution of volumetric strain. An empirical relationship proposed by Touhidi-Baghini (1998) for bitumen sands is used. This relationship reads: (2) with k 1 the updated absolute permeability, k 0 the initial absolute permeability, ε v the volumetric strain and c a constant. According to Touhidi-Baghini (1998), the values c = 5 and c = 2 appear to be appropriate to match the vertical and horizontal permeability evolutions.
The correction of the porous volume is implicitly made by using a term traducing the evolution of the porous volume, DPV, due to the geomechanical phenomena. This reads: (3) with t i the time at the beginning of the ith period, PV g k (t i +1 ) the porous volume on the FG grid deduced from the fields transferred by the geomechanical simulator at the end of the kth iteration of the ith period, PV r k (t i +1 ) the porous volume computed by the reservoir simulator and transferred on FG at the end of the kth iteration of the ith period. Furthermore, it is considered here that: where kc(i) corresponds to the iteration leading to convergence for the ith period. Furthermore, it is to note that the compressibility c p considered in the reservoir simulator is linked with the geomechanical parameters using the following relation: with b the Biot modulus, K s the matrix bulk modulus and K d the drained bulk modulus. Example of a refined grid with nested levels.

N Guy et al. / Numerical Modeling of Thermal EOR: Comprehensive Coupling of an AMR-Based Model of Thermal Fluid Flow and Geomechanics
Therefore, the Lagrangian reservoir porosity (ϕ r in Eq. 1) during the kth iteration of the ith period reads for a given gridblock: (6) with V 0 the initial volume of the considered reservoir gridblock and p 0 the initial pressure.
The Lagrangian geomechanical porosity (ϕ g in Eq. 1) reads in the present context: (7) with ε v = tr(ε) the volumetric strain related to the strain tensor ε. It is to note that the Lagrangian geomechanical porosity expression is conventional. The reservoir porosity given in Equation (5) is not conventional. The two first terms of the right member of Equation (5) correspond to a conventional expression of reservoir porosity, the third term is the correction associated with the current time period and iteration and the fourth term is the cumulated correction associated with the converged time steps. The reservoir porosity correction is set to reach the convergence criterion given by Equation (1) that ensures that porous volumes in geomechanical and reservoir simulators are consistent.

NUMERICAL STUDY
A case study aiming at evaluating the numerical algorithm is proposed herein. This case study is based on a description of Two uncoupled (i.e., considering only thermal reservoir flow) simulations are performed. The first one (RF) is achieved with a fine reservoir grid and the second one (RAMR) using AMR for the reservoir simulator. Two coupled simulations based on the presented coupling methodology are then studied. The first one (CF) is carried out with a fine reservoir grid, the second one (CAMR) using AMR for the reservoir simulation. These last two simulations are performed with the same Geomechanical Grid. In the coupled examples, CRIT is equal to 2 × 10 -3 and a period corresponds to 10 days of production. It should be also pointed out that the convergence is checked on the fine reservoir grid FG for both CF and CAMR calculations, which means that it is just as difficult to reach in both cases.

Description of the Case Study
This test case has been constructed using petrophysical and fluid properties issued from the test case presented by Zandi et al. (2010a, b). The reservoir is assumed to be homogeneous. As shown in Figure 3, the top of the reservoir is 250 m deep and the initial temperature is 10°C. The initial pressure is equal to 2.4 MPa at the top and the reservoir is assumed to be in hydrostatic equilibrium. The considered domain is rectangular with its dimensions in the X, Y and Z directions respectively equal to 150, 500 and 20 meters (Fig. 3). Here, only the half of the domain is modeled. The solution in the other half is deduced by symmetry of the theoretical problem. The dimensions of the numerical models are thus equal to 75, 500 and 20 meters. The well pair is located along the Y axis and in the middle of the theoretical reservoir along the X axis. The distance between the two wells is 6.5 m. The producer is 1.75 m above the base of the reservoir.
The simulations are performed over 1 200 days. The first 120 days consists in pre-heating the regions surrounding the well pair in order to decrease the oil viscosity and be able to inject steam and form a stream chamber afterwards in the injection well. The steam injection starts at the end of preheating with a maximal pressure set to 5 MPa. The production rate is controlled in order to keep the production well temperature 20°C to 35°C lower than the injection-well temperature. The production-well minimal pressure is set to 0.5 MPa. The steam injection temperature is about 260°C. Sideburden rocks are not modeled assuming the treated case as geometrically periodic. This assumption corresponds to the fact that in a SAGD process, it is common to use several pairs of wells that are parallel and equidistant to optimize production rates. In the model associated with the reservoir simulator, fluids cannot flow through the boundaries and heat losses, by conduction through upper and lower boundaries, is not taken into account here for sake of simplicity.
In the geomechanical model, the simulated domain includes the reservoir, surrounded below by underburden layers and above by overburden layers. The horizontal displacement of lateral boundaries is blocked as well as all the displacements of the lower boundary. Initial state stresses are supposed to be isotropic. It is to note that the considered geomechanical behaviour is taken elastic for sake of simplicity and that this hypothesis can lead to an underestimate of the contribution of geomechanical effects on production.
The mesh considered for the Geomechanical Grid (GG) is shown in Figure 4 and the domain coupled with the reservoir simulator is shown by a red rectangle. Examples of AMRbased and fine reservoir grids can be seen in Figure 5. As a comparison, when 4 000 gridblocks are used in the reservoir grid only 250 gridblocks are used to describe the coupled area in the geomechanical computation. This mesh reduction leads to a significant reduction of computational times (Guy et al., 2011) compared to the ones required by usual coupled studies using same meshes.

Results
Results obtained on the four cases as explained previously are thus compared. To characterize the global response in each case, cumulative oil produced, normalized by the total production obtained for the RF case, is plotted versus time in Figure 6.
The gap between the coupled and uncoupled simulations in terms of produced oil is of about 10%. The results clearly show that the AMR method leads to relevant results in both non-geomechanical and geomechanical contexts using the numerical algorithm proposed in the previous section. Indeed the results given by the simulation RAMR are similar to the results given by the simulation RF. The same statement can be made for simulations CAMR and CF. Figure 5 shows the oil saturation field in the reservoir after 1 200 days of production and with a coupling time step of 10 days for the two coupled cases. In all cases, one can observe that the oil saturation is low in a triangular area which corresponds to the steam chamber. The geometry of the steam chamber is the same for the two uncoupled simulations (RF and RAMR) and also very similar for the two coupled simulations (CF and CAMR). Once more the accuracy of the AMR method in both coupled and uncoupled contexts is demonstrated.
As an example of numerical results obtained with the geomechanical simulator, the vertical displacement after 1 200 days of injection is shown in Figure 7. With both simulations the top of the vertical displacement reaches 0.2 meters. Once more the results are very close with or without the AMR method. Oil saturation fields after 1 200 days of simulation for the four considered cases.
The average mesh reduction ratios between AMRG and FG for uncoupled and coupled simulations are plotted in Figure 8. The value associated with a given time corresponds to the average of the mesh reduction ratio from the start of the simulation to this time. The proposed scheme has led to a final mesh reduction ratio of about 75% considering all the simulations. As it can be seen in Figure 8, the mesh reduction is slightly lower for the coupled simulation but is still very significant.
Finally, the use of the proposed algorithm, based on AMR reservoir modeling and DAM based geomechanical coupling, instead of a classical coupling procedure based on the use of a single same static mesh, has led to a reduction of the mesh size of more than 90% for the geomechanical analysis and of more than 70% for the reservoir computations.

CONCLUSIONS
This paper presents a coupling procedure integrating both a field transfer module (based on a diffuse approximation method) and an AMR description of the thermal fluid flow. This procedure is dedicated to perform, on one hand, simulations coupling geomechanics and thermal reservoir flow with separate grids and Adaptive Mesh Refinements (AMR) to follow the flow interfaces on the other hand. The procedure has been evaluated on a synthetic but complex case describing a SAGD-based oil recovery that involves thermo-poromechanical phenomena. It has been shown that the sequential iterative coupling approach combining AMR and non-AMR contexts has led to a relevant description of the considered physical phenomena. Furthermore, the presented procedure has appeared to significantly reduce the mesh size required for the computation. The mesh reduction is of about 90% for the Geomechanical Grid (GG) in the coupled area and of more than 70% for the reservoir grid compared to classical coupled analyses without a reduction of accuracy in the computation. The proposed procedure can thus provide relevant results with significantly reduced mesh sizes.   Average mesh reduction ratios between AMRG and FG for uncoupled and coupled simulations. Vertical displacement in the reservoir after 1200 days of simulation for the two coupled simulations.