Assf Mohm-Hussin ,Plo Enriqu Vrgs Mnoz ,Polo Frnso Dloso ,Clui Sorgi ,Vinnzo D Gnnro ,Surj Kumr Suih ,Qingli Ni ,Jos´ Mri Sgur Srr ,M.R.Lkshmiknth ,Jos´ Alvrllos Iglsis ,R Dnils

a Schlumberger Geomechanics Center of Excellence,Gatwick Airport,United Kingdom

b Schlumberger,Calle Alvarez de Baena 4,Planta Baja,28006,Madrid,Spain

c Schlumberger,Av.A.Nobel,ZI Nord,Pau,64000,France

d Schlumberger,45th Floor Tornado Tower,West Bay,Majlis Al Taawon,Doha,Qatar

e Repsol Technology Center,M´ostoles,Madrid,Spain

Keywords:CHOPS Numerical simulation Sanding Geomechanics

ABSTRACT The contribution of geomechanics to provide a rigorous quantification of porosity changes and associated permeability changes is often neglected when considering hydrocarbon production.However,it often has significant effects on production rates and ultimate recovery.For the cold heavy oil production with sand (CHOPS) technique,geomechanics is the key for reservoir simulations and promoting successful operations.In fact,the technique employed in CHOPS significantly affects the stress state within the reservoir by inducing “on-purpose” formation damage.It is not only vital to comprehend the behaviour of the reservoir during production,but it is crucial to identify how to harness that behaviour to improve productivity.In order to simulate the mechanical behaviour of unconsolidated sand material,an elastoplastic damage model was formulated.Fluid-flow-geomechanical modelling was then performed for predicting individual well behaviour and overall field performance.The combined effect of fluid-flow and geomechanics improved predictions with respect to oil,water,and gas production rates at key wells.Fluid rates matched satisfactorily most of the wells.In addition,the onset and propagation of equivalent wormhole networks were quantified throughout the production history.This enabled quantification of the volume of produced sand at individual wells.The comparisons between the measured and simulated sand volume rate at well locations showed reasonable agreement.Such calibrated models can then be used for the placement of new wells to optimize production.

1.Introduction

Cold heavy oil production with sand (CHOPS) is a nonthermal technique that involves initiation of sand influx into oil wells for enhancing productivity.Modelling of CHOPS is challenging due to its complex mechanisms involving the growth of high-permeability channels (wormholes) network and foamy oil drive.It is generally believed that the growing nature of wormhole network will change the mechanical properties of formation rocks,reducing rock strength while increasing porosity,permeability etc.,therefore alter the stress conditions and flow behaviour of the reservoir dynamically[1].Also,when the pressure in the reservoir drops to bubbling point and dissolved gas is emerging from oil solution,the bubbles remain dispersed and move with oil/sand mixture as a complex multi-phase flow.This phenomenon causes the heavy oil to swell and flow with more energy.Since the free gas saturation remains zero,the mobility of oil is not impeded.The foamy oil also acts as an additional conductivity enhancing factor,leading to a very different mechanism from that of classical sand transport within porous media.

Therefore in reservoir simulation,to better understand reservoir response using CHOPS,and potentially improve the recovery factor,it is important to apply a multidisciplinary approach that integrates dynamic fluid-flow and geomechanics to enable capturing the dynamics of reservoir behaviour including the coupling mechanisms associated with wormhole network generation and enhanced fluidflow.

Fig.1.Methodology for CHOPS modelling.

Fig.2.Illustration showing non-equilibrium foamy oil and equilibrium pressure response.

The objective of this paper was to present an enhanced workflow of simulating the coupling mechanisms involved in the CHOPS process and verify modelling assumptions against production data.

In order to simulate the mechanical behaviour of unconsolidated sand material,an elastoplastic damage model was formulated.The formulated constitutive model was validated against triaxial compression tests and a sand production experiment.The numerical response was satisfactory in the prediction of mechanical features of sand material.In addition,the constitutive model was able to match reasonably the wormhole propagation reported in the sand production experiment.

In coupled fluid-flow-geomechanical simulations,the spatial and temporal variations of formation damage,shearing,and dilation were predicted by considering the onset and propagation of a wormhole network using the formulated geomechanical damage model.In parallel,the foamy oil behaviour was considered using dynamic multiphase fluid-flow simulation.In the coupling process,a finite element geomechanics simulator used the pressure gradients from a finite difference fluid-flow simulator to compute stress changes within the reservoir.These stress alterations caused deformation and damage in the reservoir rock,representing development of wormhole network with modified porosity and permeability.The porosity and permeability values were then updated in the consequent fluid-flow timestep,and the pressure was recomputed accordingly.Such method was incorporated in an automated workflow and showed the full impact of geomechanics on the reservoir behaviour and performance.

By considering rock damage and subsequent permeability enhancement,the coupled simulation improved predictions with respect to oil,water,and gas production rates at key wells in our case study,demonstrating the value of such approach in explaining and predicting well behaviour and overall field performance.Oil production rates matched satisfactorily most of the wells.In addition,the onset and propagation of equivalent wormhole networks were quantified throughout the production history.This enabled quantification of the volume of produced sand at individual wells.The comparisons between the measured and simulated sand volume rate over three years at several key wells locations showed reasonable agreement.Improved history matching raises the confidence in model’s prediction capability therefore enables us to provide better and more measured decisions for post-CHOPS field exploitation plans.

2.Methodology

In general,the CHOPS process is supported by two different but intertwined effects:

(1) Foamy oil flow

(2) Sand failure and production

Typically,sand will be produced due to shear or tensile failure caused by changes in the formation stresses due to pore pressure change.When sand fails,it generates regions of enhanced permeability,and these contribute to the increased production observed in CHOPS fields.On the other hand,the foamy oil phenomenon is a production mechanism that is responsible for reduced in-situ oil viscosity and high apparent critical gas saturation due to trapped gas bubbles in the viscous heavy oil,both of which contribute to increased production from the reservoir.

Fig.3.Foamy oil pseudo-components.

Fig.4.Illustration of dilatant and wormhole network zones around a producer well.

The methodology adopted in this paper relied on explicitly coupling the fluid-flow reservoir simulation with geomechanics simulation.Sand failure and production were modelled in the geomechanical finite element simulator where stress changes determined using the reservoir pressure changes to calculate due to the fluid-flow in the porous medium using finite difference reservoir simulator.The changes in the reservoir properties (e.g.,permeability and porosity) due to formation failure were then fed back to the fluid-flow dynamic simulator,which incorporated these effects into the fluid-flow calculations for consequent timesteps.Fig.1 shows the methodology for CHOPS modelling.

3.Chops theoretical mechanisms

As mentioned,CHOPS is controlled by two main aspects namely foamy oil flow and sand failure and production.It is thus important to formulate numerical models that can capture these mechanisms accurately.

3.1.Foamy oil behaviour

Fig.5.Concept of equivalent damaged zone.

The phenomenon called foamy oil is a production mechanism that is responsible for reduced in-situ oil viscosity and high apparent critical gas saturation,both of which contribute to increased production from the reservoir.Foamy oil is distinctive of certain heavy crude oils that exhibit non-equilibrium pressure/volume/temperature (PVT) behaviour caused by its high viscosity.In light and medium oils,when pressure falls below the bubble point,the liberated gas bubbles rapidly coalesce,and the system reached an equilibrium (saturated) condition promptly.In heavy oils,the liberated gas bubbles do not quickly coalesce;they remain trapped in the liquid (supersaturated condition) and move along with it,effectively creating a single-phase flow of oil and trapped gas.The coalescing or nucleation of gas bubbles is a diffusion-based mechanism,and in high-viscosity heavy oils,diffusion process is very slow.Compressible gas bubbles are comparatively smaller(typically smaller than pore-throat size) in foamy oil,and they do not contribute to pressure maintenance as a free gas would do in non-foamy oils.Thus,the pressure depletion of foamy oils is much steeper when compared to the equilibrium behaviour which is a typical feature of non-foamy oils.This effect is illustrated in Fig.2 and shows the pressure variation of a typical constant-rate,gas exsolution experiment response.Such response was shown by Uddin [1]in a live heavy oil sand pack conducted by Alberta Research Council.The steep pressure decline liberates more gas that eventually starts to nucleate.At this point,gas bubbles start to grow,and the high compressibility of the gas starts to increase the pressure in the system.As pressure increases,the condition of supersaturation is alleviated,slowing the driving force for nucleation.Eventually,the rate of pressure increase slows to the point where pressure starts to drop again,and the foamy oil pressure behaviour starts to approach the “equilibrium” behaviour.In foamy oils,the gas is not readily liberated and remains trapped in the oil meaning that the phase viscosity of the flowing oil is much lower.Additionally,since very little free gas is flowing,oil mobility remains high because it does not have to compete with a much more mobile gas phase.Both phenomena contribute to higher recoveries observed with foamy oil systems.

An example of damage created by gas exsolution,conceptually similar to CHOPS’s mechanical behaviour was observed by Sultan et al.[2]while investigating the mechanical behaviour of gascharged marine plastic sediments in the Gulf of Guinea.In this case,the shear strength of sediments following exsolution and subsequent consolidation is affected by the damage caused to the sediment structure by exsolution.This leads to higher compressibility due to the existence of both occluded gas bubbles in the sample and clay de-structuration.The authors highlighted that based on experimental evidence,gas exsolution and expansion induce an important degradation of the materials;to correctly simulate the experiments,a constitutive geomechanical model including damage and degradation of mechanical properties was necessary.

Regarding flow,the conventional black-oil PVT approach used in the numerical simulation is not adequate to model the behaviour of foamy oil.The model presented in this methodology employs a sophisticated treatment of fluid PVT behaviour based on kinetic reactions.This approach,developed by Alberta Research Council[1,3],uses a total of three oil phase pseudo-components (dead oil,dissolved CH4gas,and dispersed CH4bubbles) and one aqueous component(water).The exsolution of dissolved CH4is represented by the following four steps (Fig.3):

Fig.6.Evolution of cohesion and friction coefficient with equivalent plastic shear strain within the hardening and softening phases.

Fig.7.Concept of sand mobilization.

(1) Generation of dispersed bubbles

(2) Generation of connected bubbles

[gd]is dissolved CH4gas in oil phase

[bd]is dispersed CH4gas bubbles in oil phase

[bc]is connected CH4gas bubbles in gas phase (free gas)

The foamy oil behaviour is implemented into a finite difference fluid-flow simulator.

3.2.Geomechanical behaviour

Shao et al.[4]proposed an elastoplastic damage model to simulate sand production due to rock dilation and wormhole development.The proposed constitutive model used the Drucker-Prager failure surface to describe the mechanical behaviour of sands.Since the Mohr-Coulomb failure surface is most widely used in oil and gas applications,Mohamad-Hussein et al.[5]proposed an elastoplastic damage model consuming the Mohr-Coulomb failure surface,which was generally adopted and enhanced in this paper.Further modifications to the constitutive model were formulated and are explained in this section.

Sand production relied on a macroscopic erosion model.Due to stress and pressure variations,a network of wormholes started to develop and propagate around the producing wells.Fig.4 shows the dilatant zone and the wormhole network around a producer well.The heterogenous zone where wormholes were supposed to develop was considered here as an equivalent continuous damaged zone.Fig.5 shows the equivalent damaged zone concept.A scalar damage variable,d,was defined in the classic damage mechanics continua [6]and represented the reduction of effective material area due to growth of cracks:

S and Siniare the current and initial surface areas,respectively.With this definition,the damage variable was taken as the change of total average porosity,due to rock dilation and shearing during plastic flow.Because the growth process of wormholes was not explicitly determined in this approach,the onset of wormholes was supposed to be associated to the onset of plastic shear strain.Therefore,the total average porosity was expressed as:

Fig.8.The evolution of permeability change with damage.

Fig.9.Flowchart illustrating model stages.

n and niniare the current porosity and initial porosity,respectively.The term ＜εvp＞ is the dilatant plastic volumetric strain;it represents sand dilation around the wellbore.γpis the equivalent plastic shear strain;it represents the average porosity due to wormhole development.The parameter λdis a material property that describes the rate of porosity increase due to wormhole development and propagation.Accordingly,the damage criterion was expressed by:

In the material where wormholes could appear,the yield,failure,and potential surfaces were based on the Mohr-Coulomb failure criterion.The Mohr-Coulomb failure surface is written as[7,8]:

Cfand φfare the cohesion and friction angle at failure,respectively.I1′represents the first invariant of the effective stress tensor σij′,J2is the second invariant of the deviatoric stress tensor Sijand θ is the Lode angle,which describes the variation of yield stress with loading direction in the deviatoric plane.The three invariants are defined by:

The yield surface is deduced by introducing hardening law into the failure surface.The resulting yield surface is expressed as described by Mohamad-Hussein and Ni [5]:

φpand Cpare the friction angle and cohesion in the plastic hardening phase,respectively.Hyperbolic hardening laws were introduced to model the evolution of cohesion and friction coefficient with equivalent plastic shear strain [5]:

H1and H2are two hardening parameters that describe the evolution of cohesion and friction coefficient in the hardening phase,respectively.

Fig.10.Triaxial compression tests-Comparison between experimental data and numerical simulations (experimental data reported by Ref.[10]).

The elastic properties and strength properties are affected by formation damage due to sand production and expressed as [5]:

K0,G0are the initial bulk modulus and shear modulus,respectively.In contrast,Kd,Gd,φd,and Cdrepresent the bulk modulus,shear modulus,friction angle,and cohesion for damaged material,respectively.Fig.6 illustrates the variation of cohesion and friction coefficient within the hardening and softening phases.

In addition to rock failure,sand mobilization within the porous media is an important mechanism for predicting sand production.Therefore,a sand mobilization model was used in this paper.Fig.7 illustrates the sand mobilization of point A into a producer well.The sand mobilization model assumed that sands travelled with the same velocity as oil after they became loose in the damaged zone.The fluidized solid concentration was not solved explicitly,but only the macroscopic consequences of the formation of a disturbed zone were considered,i.e.,increase in average porosity and permeability,or material damage.The influence of fluidized solid on fluid-flow may be considered through variation in the sand/oil mixture viscosity.Darcy’s law was applied for the fluid mixture flow:

where k is the intrinsic permeability of the porous media,μ is the dynamic viscosity of the fluid mixture andis the pressure gradient.

The distance,D,travelled by a sand point A was computed from the sand velocity and the production time,Tf,as:

If the travelled distance was less than the distance between the sand point and perforation,then no sand production was assumed for point A.On the other hand,if the travelled distance was greater than the distance between point A and the perforation,sand volume,Vs,was produced according to the following equation:

Fig.11.Geometry and imposed pressure gradient for the sand pack.

Fig.12.Comparison between numerical simulation and digitized CT-scan experimental result (dotted red line after [11]).

Fig.13.Porosity profile along x-axis.

The increase in porosity in the damaged zone is followed by a corresponding enhancement in permeability.This permeability enhancement is used to model the role of wormholes as improved flow channels in CHOPS.In the proposed coupled approach,the new permeability calculated by the geomechanical simulator is fed back to the fluid-flow resrvoir simulator.By updating the permeability,the coupled approach allows for incorporating the effects of formation failure and sand production into fluid-flow calculations.The choice of the permeability enhancement law is an important mechanism in the coupled approach.The law should reflect the sand production behaviour.For example,the nucleation and growth of wormhole networks lead to a significant increase in the permeability.For this reason,it was decided to relate the permeability to formation damage(average porosity).As suggested by the work of Duveau et al.[9],the following third-order power function was used:

Fig.14.Overview of the 3D geomechanical model showing one horizon,reservoir pressure,and Young’s modulus.

Fig.15.Lateral resolution within the reservoir.

Where d is the damage,λ is a parameter describing the rate of increase of permeability with damage,k is the permeability and k0is the initial permeability.

Fig.8 illustrates the shape of the permeability enhancement with damage.

3.3.Flowchart of model stages

Fig.9 illustrates the flowchart of model stage.The modeling process is divided into five main stages:

(1) Stage 1:In this stage,the geomechanical model includes mechanical properties,pore pressure and in-situ stresses at preproduction condition.

Fig.16.Vertical resolution within reservoir was high to capture the vertical variation of mechanical properties:distribution of Poisson’s ratio(left)and cell thickness(right).

Fig.17.3D view showing the distribution of gamma ray.

(2) Stage 2:Fluid-flow response is obtained using a fluid-flow finite difference simulator.In this stage,pore pressure distribution within the reservoir and oil production rates at wells are computed.

(3) Stage 3:Geomechanical response is obtained as a result of the pore pressure change within the reservoir.First,poroelastic stresses and strains are computed using a finite element simulator.Then,sand production analysis is verified.If plastic deformation is induced,the damage zone is identified,and sand volumes and rates are quantified.The elastic and strength properties are degraded with damage.

(4) Stage 4:Porosity and permeability are updated as a result of sand production.Consequently,reservoir pressure and oil production rates are also affected.

(5) Stage 5:In this stage,comparisons between simulated and measured oil production rate,sand production rate and cumulative sand volume are examined.

4.Simulations of laboratory experiments

The proposed constitutive model is validated by simulating triaxial compression tests at various confining pressures and by comparing the numerical response of the model with a sand production experiment.

4.1.Triaxial compression tests

The model parameters were determined from a sequence of conventional triaxial tests performed on Athabasca McMurray oilbearing sands reported by Kosar [10].Table 1 summarizes the values of the model parameters.Numerical simulations were conducted for triaxial compression tests at various confining pressures.Fig.10 compares the numerical simulations and experimental data for the triaxial compression tests with confining pressures of 50 kPa,150 kPa,250 kPa,and 450 kPa,respectively.These comparisons show good to reasonable agreement between the model’s response and the experimental test data.The model was capable to predict the key features of oil-bearing sand behaviour such as yielding,hardening,and failure.

Table 1 Model parameters.

Table 2 Dates of selected stress steps.

4.2.Sand production experiment

Fig.18.Lithology cube in the field and within the reservoir.

Tremblay et al.[11]conducted a sand pack experiment.The sand pack had a cylindrical shape with a length of 85 cm and a diameter of 10 cm.Viscous fluid was allowed to flow through one end of the sand pack,resulting in sand production at the other end where a small orifice of 10 mm existed.The prescribed pressure at the inlet was 5.2 MPa and 2 MPa at the outlet.Fig.11 shows the geometry and imposed pressure gradient for the sand pack experiment.To visualize the formation of wormhole inside the sand pack,CT-scan pictures were taken.In this section,to show the predictive power of the constitutive model,the sand pack experiment was simulated assuming an axisymmetric problem.In the absence of triaxial compression tests for the sand pack material,mechanical properties derived from the previous section were used,and the damage parameter λdwas computed from this example.Fig.12 shows the comparison of wormhole geometry between the CT-scan experimental results and numerical simulation response,respectively.In general,there is reasonable match between the experimental and numerical results.However,the CT-scan showed a sharp contrast at the tip of the wormhole [11],whereas the numerical simulation predicted smooth transitional zone from the tip to the intact sand away from the high-porosity region as depicted in Fig.13.The predicted extent of the wormhole length was approximately 10 cm.The value of the damage parameter λdobtained from this example was 400.

Fig.19.3D property arrays describing DTC (us/m),density (kg/m3),Young’s modulus (GPa),Poisson’s ratio,friction coefficient,and UCS (MPa).

Fig.20.Average reservoir pressure (MPa) in the field between October 2014 and July 2017.

5.Field case study

The case study is devoted to modelling a CHOPS field.For this purpose,a full-field 3D mechanical earth model (MEM) was constructed.The geomechanical model considered input data collected from multidisciplinary sources,such as seismic data,well logs,and fluid data.The model was then calibrated in detail with field measurements to match oil and sand production data in individual wells at various stages of reservoir depletion.A reasonably history matched model is expected to significantly increase the predictive power of the model and provide guidelines for the placement of new wells in regions enabling production enhancement.

Fig.21.Pore pressure distribution (MPa) within the reservoir at four-time steps.

Fig.22.Total vertical stress distribution (MPa) within the reservoir at four-time steps.

Fig.23.Total minimum horizontal stress distribution (MPa) within the reservoir at four-time steps.

Fig.24.Distribution of formation damage within the reservoir at four-time steps.

Fig.25.Wormhole development and propagation within the reservoir zone.

5.1.3D mechanical earth model

Data from 245 wells were used for horizons construction,property modeling,and validation.There is excellent coverage of well logs,which included density,gamma ray,porosity,and compressional slowness along these wells.Among the 245 wells,production data such as oil,water,and gas and sand rates and volumes were measured for 8 wells situated in the zone of interest.In addition,37 wells had formation tops collected from geological and drilling reports.There were 3 wells with a reasonable set of well logs that contributed to 3D MEM construction.There was not enough data to fully calibrate the maximum horizontal stress magnitude in the study area.However,there was relevant maximum horizontal stress magnitude reported by Bell et al.[12]within the region of the studied field.The interpretation was based on two micro-fracturing tests conducted at two depths in the Deadwood formation in Regina (well 3-8-17-9W2).The average maximum horizontal stress and minimum horizontal stress ratio(σHmax/σhmin) that was computed at these two intervals was 1.33.

A 3D Mechanical Earth Model was constructed for an area of approximately 100 km2(Fig.14).The 3D grid was built for the entire reservoir.The lateral resolution was chosen to be 20m by 20m around the eight key wells with production data(Fig.15).Twentyfive horizons were used to define the geometries of the overburden and underburden layers.The vertical resolution of the model was approximately 30 cm which is considered very high and able to capture the detailed vertical property variation (Fig.16).The 3D grid was populated with mechanical properties derived from well logs,and the preproduction pore pressure was from the reservoir fluid-flow model.The first step for populating the mechanical properties was to define a 3D facies cube.Gamma ray was the basis for deriving the lithology cube:

(1) Gamma ray logs were available on most of the wells.

(2) Normalized Gamma ray logs were upscaled from the well logs into the 3D grid.

(3) The upscaled gamma ray was populated geostatistically between the wells,as shown in Fig.17.

(4) The following gamma ray cut-offs were used to distinguish sand from shale lithologies:-100 API for over burden zones.-80 API within the reservoir

Fig.18 shows the obtained lithology cubes within the reservoir and surrounding formations.

For property modelling,the mechanical properties such as Young’s modulus,Poisson’s ratio,unconfined compressive strength (UCS),friction angle and tensile strength were upscaled from the 1D MEMs into the 3D geomechanical grids along the well.The upscaled properties were then populated using the compressional slowness (DTC) cube as co-kriging.The following correlation was used to convert dynamic to static Young’s modulus following the relationship Esta=0.032Ed1.yn632.The static Poisson’s ratio was assumed to be equal to the dynamic Poisson’s ratio.Fig.19 shows the 3D cubes of the individual mechanical properties.

Fig.26.Distribution of the permeability ratio (k/k0) within the reservoir at four-time steps.

Fig.27.Distribution of pore volume increase (m3) within the reservoir at four-time steps.

It is worth to mention that the wormhole network occurs only in the reservoir where fluid flow is assumed active.Its development and propagation weaken formation rock and induce deformation at various locations within the reservoir.The numerical model computes dynamically the changes in porosity and thus updates the mechanical properties based on such changes.

Fig.28.Comparison between reservoir simulation results of the coupled (solid line),uncoupled simulations (dotted line),and the observed (circles) field oil production rate.

The movement of the horizons is governed by the magnitude of deformations obtained by the wormholes.Such movement exists around the wormhole network.The dynamic movement of the horizons is captured in the finite element model.

A preproduction stress state was obtained considering the gravity and pore pressure forces acting internally on the model and an external horizontal loading scheme acting on the model boundaries to represent the maximum and minimum horizontal stresses.Displacement boundary conditions were used for the external horizontal loadings.The advantage of this approach is that a complex stress state can be calculated within the model representing the initial effective stress state in a realistic manner because it accounts for the geomechanical property distributions and the shape of the reservoir.A strike slip dominated stress regime was adopted in this case study as reported by Bell et al.[12]for the Western Canada sedimentary basin.In addition,a maximum horizontal stress direction of N50.9°E was used as reported by Hamid[13].The preproduction stress state was used as an initial stress state for the production modelling phase.During production,these stresses will change,and their effect on wellbore behaviour and sand production must be reassessed.Moreover,the changes in the stress might cause rock deformation,which will modify its porosity and permeability.This,might change the fluid production and sand production rates at the wells.

The production scenario in the dynamic fluid-flow model started in October 2014 and ended in July 2017.A stress step scheme was required to determine the optimum timing to incorporate the effects of oil/gas/water production on sand production.The stress steps were the points in time at which stress analyses were performed,accounting for maximum fluid effects,and providing suitable points in time for verification of geomechanical-related events such as sand production rate.In this paper,12 stress steps were selected,as summarized in Table 2.

5.2.Results

Fig.20 illustrates the variation of average reservoir pressure in the field over time.The curve shows a decline and increase in pore pressure.The decline in reservoir pressure liberated more gas that eventually started to nucleate.At this point,gas bubbles started to grow,and the high compressibility of the gas started to increase the pressure in the system.Fig.21 shows the pore pressure distribution in a representative reservoir layer at different stress steps.The purple zones localized around the individual producer wells represented the regions of high production.In the coupled simulations,the changes in pore pressure are due to productioninduced changes in effective stresses.The greatest changes in effective stresses were in regions around the producing wells.Fig.22 and Fig.23 show the distribution of total vertical stress and total minimum horizontal stress in the reservoir,respectively.The total vertical stress was controlled by the densities of different lithology and thus the coupled modelling results showed slight changes in the total vertical stress.On the other hand,oil production resulted in significant changes in the total minimum horizontal stress.If the changes of the total and effective stresses are significant,it could cause formation damage that might lead to sand production.The distribution of formation damage is presented in Fig.24.Formation damage represents the average porosity increase due to development and progressive propagation of wormholes around producer wells.The largest damage zones occurred at well locations and the magnitude declined away from the producer zones.The development and propagation of wormhole networks is illustrated along a vertical cross section(Fig.25).Note that the reservoir is a porous medium and,in modelling the coupling between geomechanics and fluid-flow,various factors govern the fluid flow.The permeability is a principal factor that governs transport of fluid (oil/water/gas) to the surface via the wells.The wormhole network development and progression resulted in permeability enhancement.The distribution of permeability ratio k/k0at various times within the reservoir is presented in Fig.26.The pore volume increase due to sand production is shown in Fig.27.These figures illustrate that the permeability increases due to sand production and formation of wormhole networks within the reservoir.This permeability change lead to a pressure gradient by which the oil flow in individual wells increased substantially.It is important to note that the general trend of flow and the development of wormholes were aligned along the regional maximum horizontal stress direction.Although further study is required to establish such a connection,this association provided additional evidence that geomechanics was an inherent component of reservoir behaviour.This permeability change led to a pressure gradient by which the oil flow in individual wells increased substantially.

It is important to note that the bottom hole pressure measurements were not available.Instead,the tubing head pressure were recorded at the wells.The tubing head pressure was extrapolated to the reservoir level to compute indication of the bottom hole pressure at the study wells.The derived bottom hole pressure was then used as a guidance in the fluid flow simulations.Hence,the emphasis was to achieve history match of oil and sand production rates.

Oil production rate,sand volume production rate,and cumulative sand volume production were plotted for key wells.The comparison between coupled model prediction and actual oil production rate in the field is presented in Fig.28.The comparisons between simulated and observed oil production rate in the key wells are shown in Fig.29.These comparisons show overall improvement in well performance compared to original fluid-flow simulation.The impact of geomechanics was clearly demonstrated throughout these comparisons.An important component in the proposed approach to model CHOPS was to properly simulate sand production.To predict sand production,a sand mobilization algorithm implemented in the geomechanical simulator was activated.The plots shown in Fig.30 and Fig.31 depict the total cumulative sand volume and sand volume rate in the wells,respectively.These comparisons show good agreement between simulated and historical sand production.The model was able to predict high sand production rate and low production rate at the various wells.

Fig.29.Comparison between reservoir simulation results of the coupled(solid line),uncoupled simulations(dotted line),and the observed(circles)oil production rate in key wells.

5.3.Guidelines for well placement

The results of the 3D geomechanical simulation can provide valuable guidelines for the placement of new wells.Fig.32 shows the distribution of formation damage and permeability ratio within a representative reservoir layer for the last considered time step.The following general guidelines were deduced,which naturally should be considered within the further field development strategy:

(1) Oil producer wells (OP wells) are highlighted by black ellipses.Placement of new wells within zone of high permeability will increase the chance of producing high volumes of oil compared to sands.

(2) CHOPS wells are highlighted by white ellipses.Placement of new CHOPS wells is generally recommended in zones with low permeability and wormholes.

6.Conclusions

An elastoplastic damage model was used to predict sand production associated with cold heavy oil production.The elastoplastic damage model was based on the Mohr-Coulomb failure surface.Sand production was based on equivalent damaged zone concept that quantified rock dilation and wormhole development and propagation.The damage variable was taken as the change of total average porosity,due to rock dilation and shearing during plastic flow.

Fig.30.Comparison between predicted (solid line) and observed (circles) cumulative sand production in key wells.

A coupled fluid-flow-geomechanical model was constructed for a CHOPS field with the following workflow:

(1) The model grid was constructed for the reservoir and surrounding formations.Population of the 3D mechanical properties was based on the 1D MEM properties.

(2) The reservoir behaviour was governed by the elastoplastic damage model whereas the surrounding formations followed an elastoplastic behaviour.

(3) By connecting the two independent simulations for fluid flow and geomechanics,the process of history matching was performed with a well-defined permeability updating model.

(4) The permeability enhancement law used in the model included permeability/damage relationship that controlled the effective permeability of the equivalent wormholedamaged zone.

(5) The non-equilibrium behaviour exhibited by foamy oil considered the combined effects of sand failure,mobilization,and production.

Fig.31.Comparison between predicted (solid line) and observed (circles) sand production rate in key wells.

The main results of the study show:

(1) The coupled simulation showed good match in measured oil production rate at field level.This is a marked improvement over the results of non-coupled simulation where geomechanics was not considered.

(2) History matching of oil production rates at six key wells also showed similar improvement over the non-coupled simulation.In fact,the simulation results matched the real well data at most wells.

(3) Both accumulated sand production volume and sand production rate also closely matched measurements at six key wells.

(4) It is important to note that the general trend of flow and the development of wormholes were aligned along the regional maximum horizontal stress direction.

(5) Based on the permeability change,it is recommended that future oil producing wells should be in the area with higher permeability,while CHOPS wells should be in the areas with low permeability.Both should avoid existing extensive wormhole network.

This study clearly demonstrates that geomechanics is an inherent component of reservoir behaviour.Using suitable damage and permeability enhancement model to capture rock behaviour more accurately,with proper calibration of coupled fluid-flowgeomechanical simulation is a viable tool in understanding reservoir behaviour and providing post-CHOPS production plan such as placement of new wells.

Fig.32.Distribution of formation damage (left) and permeability ratio k/k0 (right) at present day.

Declaration of competing interests

No conflicts of interest associated with this publication.As corresponding author,I confirm that the manuscript had been read and approved for submission by the authors.

Acknowledgments

The authors thank Repsol and Schlumberger for permission to publish this paper.