Particle Scale Reservoir Mechanics

Particle Scale Reservoir Mechanics — This paper summarizes some developments and applications of discrete particle modelling based on the commercial code PFC (Particle Flow Code). The following examples with relevance to geomechanical reservoir behaviour are described: – Numerical and experimental results are shown illustrating nonlinear evolution of depletion-induced compaction and stress path for simulated reservoir rock, i.e. materials that are formed under stress. This nonlinearity is often lost in core measurements, as a result of stress release during coring (core damage). – Angular and breakable superparticles have been implemented in the particle model and applied to simulate the formation of a compaction band at high stress level. The model is also used to investigate the formation of a shear band at lower confinement. – Numerical scheme of fluid coupling has been developed in both 2D and 3D models and applied to study stress-dependent permeability. Oil & Gas Science and Technology – Rev. IFP, Vol. 57 (2002), No. 5


INTRODUCTION
Reservoir geomechanics has become an area of increasing focus within petroleum reservoir management.Reservoir mechanics deals with: -Reservoir rock deformability, which may cause compaction and associated possible well problems (like casing collapse) or environmental problems (like surface subsidence; e.g.Geertsma, 1973;Chilingarian et al., 1995).
Compaction is on the other hand a drive mechanism, and may improve recovery.-Reservoir stress path, i.e. how rock stresses change as a result of production associated processes like depletion and enhanced recovery operations (Teufel et al., 1991;Addis, 1997;Morita et al., 1989).The stress path controls to what extent the reservoir deforms elastically, or if some kind of rock failure may occur during production.-Reservoir flow performance, i.e. how stresses and stress alterations affect preferred flow directions (Heffer and Dowokpor, 1990) and reservoir permeability (e.g.Holt, 1990;Ruistuen and Hanssen, 1996;Boutéca et al., 2000), and the coupling of fluid flow simulation to geomechanical models (e.g.Gutierrez and Lewis, 1998).The failure mechanism and the constitutive behaviour (in particular dilatancy vs. contractancy) have a profound impact on permeability.-Reservoir monitoring, i.e. the use of repeated seismic surveys ("4D seismics") or the use of permanently deployed geophones to monitor production-induced changes in the reservoir.
In this paper, we will discuss all these aspects of reservoir geomechanics, but from the very small-scale point of view.We will demonstrate how the use of a numerical tool, namely discrete particle modelling, can provide insight into particle scale mechanisms that may be important for the understanding of complex reservoirs.We use PFC (Particle Flow Code) as our platform for these simulations.Our numerical work is strongly correlated with laboratory experiments, which are presented in other publications.Furthermore, we do not attempt to upscale our results, but will refer where appropriate to related field-scale observations.
PFC is a simplified implementation of the distinct element method (Cundall and Strack, 1979;Itasca, 1999a and1999b).It simulates an aggregate of cemented or uncemented disks (2D) or balls (3D).The particles are rigid but can overlap or detach.The contact force between two contacted particles is determined from the overlap and relative movements of the particle pair according to a specified contact law.The particle motion is governed by Newton's second law.Particles may be bonded together in order to simulate e.g.cemented sandstone.The bond breakages mimic micro damages developed in the sandstone samples.Designed with the discrete feature, such a model can properly simulate a granular material with inherently discontinuous properties.

RESERVOIR ROCK DEFORMABILITY AND STRESS PATH
A main difference between rocks within a reservoir and on the Earth's surface is the difference in stress and stress history.While we can apply in situ stresses to a reservoir rock core in the laboratory, we can usually not compensate for possible effects of the stress history.An example of this is core damage, induced by the stress release occurring as the core is drilled.Holt et al. (2000a) presented laboratory simulations where this problem was systematically studied by synthetic rocks manufactured under applied stress.Holt et al. (2000b) showed that the main features of these laboratory experiments could be reproduced by use of PFC.
Figure 1 shows a typical stress-strain curve for uniaxial (oedometric) compaction of a granular material cemented under stress.This is a PFC 3D simulation with a granular assembly of 31 617 spherical particles, with a prescribed size distribution resulting in 20% porosity.Parallel bonds were inserted in the sample with initial axial stress 30 MPa and lateral stress (isotropic) 15 MPa.Input parameters are fitted to give good agreement with experimental simulations (Holt et al., 2000a) with analogous synthetic sandstone formed under the same stress conditions.The shear and tensile bond strengths were chosen equal to 8 MPa, with a standard deviation of ± 8 MPa.The uncemented bond normal stiffness was chosen as 4 GPa, and the additional contribution of the bonding was also set to 4 GPa.The ratio between normal and shear stiffness was set equal to 1 for the uncemented bonds, and reduced to 0.5 as a result of cementation.
A typical feature of the observed stress-strain curve is nonlinearity, with a high initial stiffness.The material softens as it compacts, as a result of an increasing amount of intergranular bond breakage.From cyclic loading (not shown in the figure) it is found that the initial behaviour is largely elastic, whereas the material at high stresses is more plastic.The cumulative number of broken bonds (also shown in Figure 1) indicates the emerging plasticity.As can be seen, there is not a well-defined onset, but the material appears to have a plastic strain component all the way from its virgin -Stress-dependency of the acoustic properties of a bonded particle assemblage has been studied.
The preliminary results of the simulations with various developments show qualitative and also quantitative consistence with experiments.It proves the feasibility of applying a discrete particle model as a numerical laboratory to study particle scale reservoir mechanics.

Figure 1
Axial stress and cumulative number of broken bonds vs. axial strain in a PFC 3D simulation of uniaxial compaction, starting from the state of stress (axial: 30 MPa; lateral: 15 MPa) where the material was formed.
state.This is in agreement with observations in laboratory tests with synthetic sandstones (Unander and Holt, 1998).
The elastoplasticity of rocks is one of the fundamental reasons for observed difference between static and dynamic moduli.
Figure 2 shows the lateral stress development (the stress path) during the numerical uniaxial compaction experiment above.There is a relatively slow increase in lateral stress with compaction during the initial part of loading, while the lateral stress increases much more rapidly at high stress.This pattern is seen in all numerical and experimental simulations performed.The quantitative behaviour of course depends on the details: in laboratory experiments, the degree of nonlinearity is much larger in weak than in competent sandstones.In the numerical simulations, the lateral stress response increases with increasing ratio between normal and Figure 2 Lateral stress vs. axial stress in a PFC 3D simulation of uniaxial compaction (cf.Fig. 1), starting from the state of stress (axial: 30 MPa; lateral: 15 MPa) where the material was formed.
shear contact stiffness, with inter-granular friction, and with change from spherical particles to non-spherical superparticles (ref.Section 2).
Field observations are rarely done with sufficient accuracy (Holt, 1999) to confirm or disconfirm the existence of nonlinear stress paths, or nonlinear stress-strain behaviour.There are several examples that subsidence tends to accelerate during depletion.This may be caused by enhanced compaction (as predicted by particle scale modelling), but also by time-delayed consolidation of low permeability (e.g.shale) layers within the reservoir, or by stress redistribution in the overburden.

Approaching Grain Angularity
In many cases, the real grains are nonspherical.PFC approximates them with circular or spherical particle elements.The simulation qualitatively agrees with observations of complicated behaviours of granular material with the simplification.However, the ability of reproducing experimental results by particle modelling is reduced.Particle shape affects particle behaviour of interlocking and rotation.It will affect the strength, in particular the frictional component of it, in addition to deformation and porosity of the particle packing.This can be understood by intuition, and thorough studies have been published (Jensen et al., 1999(Jensen et al., , 2001)).
To properly simulate angularity of some granular material, elliptical (Ting et al., 1995;Ruistuen, 1997) and superquadric (Rege, 1996;Preece et al., 1999) particle models have been studied.If applying more general distinct element method (Cundall, 1971(Cundall, , 1988;;Hart et al., 1988) or discontinuous deformation analysis (Shi and Goodman, 1985;Shi, 2001), polygonal particles can also be implemented.The price of creating complicated particles is computational time and memory.It should be considered with the balance of accuracy.

Cluster and Clump
Angularity can be approached with a method different from those mentioned above.Super particles with different shapes can be created by several circular (2D) or spherical (3D) particles.In PFC, they are called cluster or clump.A cluster consists of several particles bonded by intragranular bonds.It behaves as one particle and is deformable.The deformation properties are controlled by the contact law and the properties of intragranular bonds.A cluster can also be bonded to neighbouring clusters with intergranular bonds.To simulate real granular material, certainly, different bond properties (stiffness, strength, geometry) should be assigned to intergranular and intragranular bonds.The particle cluster can simulate grain crushing because the intragranular bonds can be broken.A clump is slightly different from a cluster.In a clump, the intragranular bond can be thought as infinitely stiff and strong.Actually, the calculation for the contacts inside the clump is skipped, and just relatively fixing the particles inside the clump.The clump approach is more efficient than the cluster approach.The clump is rigid with deformable external contacts.It can also be broken to simulate grain crushing if introducing a breaking criterion (Li and Holt, 2001a).The advantages of the cluster approach are the deformability of the superparticle and the direct way to simulate grain crushing.The superparticle model in this paper applies the cluster approach.Figure 3 shows what a super particle packing looks like.The shapes of the superparticles in the model are irregular but restricted with aspect ratio (longer dimension/shorter dimension) < 1.5.The size distribution of the particles can be controlled.In this model, the size distribution is Gaussian in φ scale, where φ is defined as -ln d/ln2 if the grain size d is in millimetre.
Figure 3 Packing of superparticles, with irregular shapes and Gaussian size distribution in φ scale.

Shear Band and Compaction Band
Biaxial compressive tests have been done with the 2D superparticle packing under different constant confining stresses.With low confining stress, the sample fails with a band formation that is in a small angle (< 45º) with maximum principal stress.In Figure 4, the confining stress is 4 MPa.It is a low confining stress relative to the uniaxial compressive strength that is 9.9 MPa for this sample.The angle between the failure band and maximum principal stress is about 35º.The band consists of mainly tensile intergranular bond failures and a few intragranular bond failures (grain crushing).With increasing confining stress, the angle between the failure bands and maximum principal stress tends to increase.In Figure 5, the confining stress is 18 MPa.The failure band is nearly normal to the maximum stress.Shear intergranular bond failures and grain crushing are much more than in the failure band in previous case.We think that in the case of low confining stress the failure bands are shear bands, in the case of high confining stress the failure bands are compaction bands.The compaction band is defined as a failure zone normal to maximum principal stress.It may occur in porous sandstone in some cases.Grain crushing and pore collapse happen in the zone and induce critical permeability reduction.The observations of the numerical simulation are consistent with those observed in laboratory experiments (Olsson, 1999;Olsson and Holcomb, 2000).
Considering a sample of porous material with a lateral confining stress σ 3 applied, the maximum vertical stress sustained by the sample, according to Mohr-Coulomb's criterion, is: (1) where φ is friction angle and UCS is uniaxial compressive strength.This equation is valid if the grains in the sample can sustain the contact forces.The sample fails if the differential stress is large enough to satisfy Equation (1).The failure will probably form a shear band, with mainly tensile bond failures and grain rotations in the granular sample.If the maximum principal stress goes up with confining stress increase, the contact forces finally break the grains at some contacts and cause grain crushing.Since there is plenty of free space for the yielding of the debris of the crushed grains in the porous material, the crushed grains suddenly lose their support capacities.This will induce contact forces increase in the vicinity and induce more grain crushing.The grain crushing propagates along the plane normal to the maximum principal stress.If this explanation is correct, compaction bands could occur under the conditions as below: -The porosity of the material is high enough so that there is enough free space for the yielding of the crushed grains.This is also a cause of inhomogeneous stress distribution that is favourable for the onset of grain crushing.
Figure 6 Failure envelope (determined at the critical points where the sample turns to dilate or compact) of a numerical sample (Li and Holt, 2001a).
-The confining stress is high enough so that the sample will not fail in shear mode due to differential stress.-An end cap exists in the failure curve of the material (Fig. 6).The shear band forms in dilatant mode corresponding to the left part and the compaction band forms in compactive mode to the right part.

Basics of Fluid Coupling in Particle Model
Many researchers have studied fluid coupling in 2D particle models for different purposes with different approaches (Thallak et al., 1991;Bruno, 1994;Cundall, 1999;Cerasi et al., 2000;Bruno et al., 2001;Li and Holt, 2001b).Based on the authors' previous study, the algorithm of fluid coupling in the 2D model has been extended to the 3D model.A fluid coupled particle model is a combination of a discrete particle model and a flow network model.The pore network is constructed with the pore space in the particle assemblage, and will be updated with the deformation and damage developments in the solid frame.The impact of the liquid phase on the solid frame will also be counted by applying the pore pressure on the particles.The flow in the pore network is solved with a time-stepping approach.The pore system in the particle model consists of pores connected by pipes.
Here, the wide spaces in the pore system are called pores.
The throats between the pores are called pipes 1 .

Construction of Flow Network
Both in 2D and 3D, the network models are constructed in two steps.In the 2D model, first, every 3 neighbouring particles create a pore enclosed by a triangle (Fig. 7a).Each pore has 3 pipes corresponding to the 3 sides.The pore volume is related to the area enclosed in the 3 particles.
Limited to the 2D model, the flow path is not visible.A virtual flow path is assumed as in Figure 8. Assuming there is a third virtual particle outside the 2D model plane and neighbouring with the two particles in the model, the three particles form a pipe.Such a pipe is not realistic.However, it is a reasonable approximation of the real granular material in 3D space.Then, some of the pores will be merged into pore clusters (enclosed by polygons, Fig. 7b).The condition for merging two pores is that the pipe conductance is larger than a predefined value.The pipe conductance will be defined in Section 3.4.The pore mergence will improve the computation efficiency.The pore number for calculation is reduced.The widely open pores are avoided, or else it will easily take in (1) In some literatures on network model, they are named as sites or cells and bonds respectively.However, in particle model the term "bond" has been reserved for cements between particle pairs.2) has to be used for the calculation.The mergence also reduces the accuracy of the calculation.The error is related to the conductance threshold to make mergence.
In the 3D model, every 4 neighbouring particles are found first, and a pore is created between them.The pore is enclosed in a tetrahedron (Fig. 9a).The model space is meshed by tetrahedrons (Fig. 9b).The pore volume is related to the space between the 4 particles.Each pore has 4 pipes corresponding to the 4 faces.Then, similar to the scheme in 2D, some of the pores are merged.Without the assumption of the virtual pipe, the flow network model is directly constructed in 3D model.
The flow resistance of the pipe is mainly controlled by the minimum section of the pipe (Bryant et al., 1993).This section, therefore, is mainly concerned in the calculation.As in Figure 8, the section area is: (2) where S triangle is the area of the triangle, S solid is the area occupied by the solid phase.It may include the fan areas of the particles and the cement part if it should be counted.Both the pore volume and the section area of the flow path can be reduced by overgrowth of the cement.This makes the model able to cope with material with any porosity.

Update of Flow Network
The scheme of constructing the flow network in two steps is flexible.It is convenient to update the flow network in coupling calculation, with the deformation and damage developments in the particle model.If the triangle mesh is limitedly distorted with the deformation of the particle assemblage, the data structure of the pores does not need updating.Only the parameters of the pore system need to be recalculated.This can be applied in the study of stressinduced permeability alteration that will be presented in this paper.However, the pore clusters should be remerged according to recalculated pipe conductances that are changed with the deformation.There is not much computation for this task.

Pipe Conductance
The micromechanism of the pore flow is assumed to be consistent with Hagen-Poiseuille's pipe flow.The conductance of a pipe can be written as: (3) where α is a conductance factor, assumed as a constant in the calculation.It is related to grain shape, pore shape and mineral constituents, etc.It can be adjusted by comparing with experiments.r e is related to the shape and the area of the section, as an effective value.However, to be efficient in calculation, it is simply calculated from an equivalent radius of the section area.r av is the average radius of the particles forming the pipe, instead of the pipe length in Hagen-Poiseuille's Equation.The simplification with r e and r av induces error.It is compensated in the calibration of conductance factor α.
The effects of the bond breakages and grain crushing are treated as the reduction of r e by increasing the solid area in Equation (2).

Flow Calculation
The network is solved in an explicit time-stepping scheme (Li and Holt, 2001b).The fluid exchange through a pipe in a time step is calculated according to: (4) where ∆p is the pore pressure difference between two pores connected by the pipe, ∆t is the time step for fluid calculation.
The pore pressure in every pore is updated according to: (5) where p 0 and p 1 are previous and updated pore pressures of the pore, K f is the fluid bulk modulus, ∆V f is the total fluid taken in the pore in a time step by summing up ∆V f in Equation ( 4) through all pipes connected to the pore, v pore is the pore volume.Note that in this paper, the compressive deformation of the particles due to fluid pressure is not considered.This is not consistent with Biot's theory.This limitation may be avoided by introducing a scheme to update the solid compressive deformation.In that case, Equation ( 5) is still in the same form but with solid deformation included in v pore .The time step for the fluid calculation should be small enough.A critical time step is: (6)

Virtual Pore Volume for Quasi Steady State Flow
In some fluid coupled problems, for example, in the study of stress-induced permeability alteration, the pore pressure distribution is relatively stable.It can be approximated by steady state flow.The calculation efficiency can be improved by introducing a virtual pore volume concept.When the flow is in steady state, the pore pressure distribution and flow rate are independent of the pore volume.They only depend on the conductances of the pipes.If the volume of a pore is small and it is well connected to its neighbour, a small time step should be used concerning this pore because the pore pressure is easy to change.If the volume of a pore is large, on the contrary, and it is poorly connected to its neighbour, a large time step should be used.The actual calculation should use the minimum time step according to Equation (6).Then the poorly connected large pore needs more cycles to reach proper pore pressure so that the steady state flow can be reached.Instead of using the real pore volume, virtual pore volume can be assigned to a pore as: (7) Then a well connected pore is assigned a large pore volume.A poor connected pore is assigned a small volume.Therefore, all the pores need similar time to respond to pore pressure change.It saves the time for calculating the real pore volumes and makes the pore pressure reach steady state quickly.

Background
Permeability estimation is important in the management of hydrocarbon reservoirs.Reservoir stresses change during production.Many researchers have found that permeability depends on stresses (Zoback and Byerlee, 1975;Fredrich, et al., 1993;Bruno, 1994).In some cases, permeability may be dramatically changed due to damage developments (especially some critical cases like shear or compaction localization, Holt, 1990;Ruistuen, 1997;Boutéca et al., 2000).
To predict the effect of stress change on permeability, Seeburger and Nur (1984) developed a network model with the consideration of a hydrostatic stress state present.The geometric parameters of the flow network are related to the hydrostatic stress.Bruno (1994) implemented a network model in a discrete particle model to study stress-induced permeability anisotropy in sedimentary rock.The permeability alteration is a result of the particle rearrangements in the model.Even without addressing the issue of hydromechanical coupling, some of the efforts to reconstruct the network model are close to the fluid coupled particle model presented in this paper.Bryant et al. (1993) constructed a network model based on Finney's packing (1970).Finney's packing is a porous medium consisting of 25 000 balls with equal size.Bryant used Finney's measurements of the spatial coordinates of 3367 sphere centres and constructed network models.By adjusting the coordinates and the radii of the spheres, compaction and overgrowth cementation have been simulated with permeability alteration.
The fluid coupling scheme presented in previous sections is developed for studying stress-dependent permeability of a granular material.Some tests have been done.They verify the fluid coupled particle model to some extent.The Kozeny-Carman equation is used as a reference since it is often a successful estimation of permeability in porous material.The Kozeny-Carman equation can be written in the form: (8) where d is the size, ϕ is porosity, γ is a pore shape factor and T is tortuosity.An empirical estimation of γT is around 5.

Setup for Studying Stress-Dependent Permeability
Figure 10 shows the setup for the study of permeability alteration under stress in 2D and 3D model.The pore pressures on the top and on the bottom are defined and kept constant, but with a constant difference between them.The difference of the boundary pore pressures drives the fluid flow through the sample.The fluid exchanges through the lateral boundaries are assumed zero.
The model simulates monophasic flow.Before loading, the sample is assumed to be saturated with fluid, and the model is solved until the pore pressure distribution is in a steady state (Fig. 11).The pore pressure is linearly distributed along the flow direction.This is consistent with Darcy's law.
Because the variation of the pore pressure distribution in the sample is small during loading, an approximation is applied in flow calculation.The flow network is solved to  steady state, and then updates the effects of the fluid pressures on the particles and calculates the sample permeability after every n particle motion calculations.The macroscopic permeability of the whole sample is calculated according to Darcy's law, and normalized relatively to the initial permeability before loading.

Permeability Alteration Induced by Deformation under Hydrostatic Loading in 3D Model
A 3D fluid coupled particle model under hydrostatic loading is studied.The model consists of 3185 uncemented particles with diameters around 1.1 mm.The constructed flow network consists of 17 763 pores and 34 505 pipes.The porosity of such a packing is 39.3%.The initial permeability of the packing quantitatively fits the Kozeny-Carman's Equation if the conductance factor α (Eq.( 3)) in numerical model is chosen as 0.82.When the sample is compressed under hydrostatic loading, the permeability reduces.Since no damage effect is considered in the model, the permeability alteration is the effect of a hydrostatic deformation and porosity change.The permeability reductions predicted by PFC and by Kozeny-Carman's Equation are quite consistent (Fig. 12).Since the initial permeability has been adjusted by the conductance factor, there is a quantitatively good fit.

Permeability Alteration and Anisotropy in 2D Model
It is much more efficient to work with a 2D model than with a 3D model.However, 2D discrete element methods, including distinct element method and discontinuous deformation analysis etc., are not able to accurately simulate 3D problem.The main difficulties in particle model are how to exactly represent porosity and material constitutive relationship.On the issue concerning stress-induced permeability alteration, current 2D model is not recommended for quantitative analysis.
Even though the representability of porosity in a 2D model is doubted, Figure 13 shows a good consistence between Kozeny-Carman's Equation and PFC prediction, if the 2D model can be assumed as a slice of granular material under plane strain state and the initial porosity can be assumed to be 30%.The 2D model was originally designed to simulate a real glass beads packing with porosity 36-37%.The diameters of the glass beads are 1.00-1.18mm.In the figure, porosity is calculated according to volumetric strain by assuming the 2D model a plane strain state and the initial 533 Figure 13 Permeability alteration under plane strain state with loading equally in x and y directions.Comparison of PFC 2D and Kozeny-Carman's Equation.The initial porosity of the sample is assumed 30%.The porosity alteration is calculated from the initial porosity and volumetric strain.The 2D model is applied to study stress-induced permeability anisotropy (Fig. 14).The sample is loaded in directions normal to or along the flow direction.When loaded normal to the flow direction, the permeability of the sample decreases with stress increase.After the sample fails, the permeability rebounds a little due to the dilatant deformation.When loaded in the flow direction, the permeability of the sample decreases very little at the beginning.Then the permeability keeps constant or slightly increases until the failure point.After the sample fails, the permeability increases a lot due to dilatancy.The calculation by the Kozeny-Carman's Equation according to porosity change in the sample does not distinguish the difference.
The effect of grain crushing has been mentioned in Li and Holt (2001b), where the critical permeability reduction is simulated (Fig. 15).The effect of different bond breakage modes (shearing or tensile) has been studied.We expect to publish some results on that issue soon in a separate paper.

534
Figure 15 Dramatically compactive deformation (a) and permeability reduction (b).K is stress path defined as the ratio of increment of vertical stress and lateral stress (Li and Holt, 2001b).

Basics of the Simulation
PFC is a dynamic model.Wave propagation can be directly simulated.From a packing of particles (as in Fig. 16), a few particles are chosen as a transmitter and some particles or groups of particles away from the source are chosen as a receiver or a receiver array.Several receivers are useful to study the wave propagation in more detail, and to get a reliable estimate of velocity.Also, the receivers need to be sufficiently far from the transmitter that near field effects can be neglected.To simulate the wave propagation, a sine velocity pulse is applied to the source particles.
Figure 16 In a particle model, some particles are chosen as transmitter and some as receivers.
The disturbance of the particle velocity will propagate in the particle packing (Fig. 17), either as a P-wave or as an S-wave, depending on the excitation of the transmitter particles.The average velocity of the particles in a receiver is monitored.With the wave travel time and travel path, the phase velocity of the wave can be calculated, in exactly the same way as in laboratory tests.It was checked that the "dynamic" elastic moduli calculated from the measured velocities were identical to "static" moduli measured directly in small stress cycles within the elastic range.A benefit of the numerical method is that the transmitter and receiver can be placed inside the sample, which reduces problems of coupling the elastic energy and surface reflections.

Examples of Velocity Measurements with PFC
Figure 18 shows examples of wave velocities recorded during simulations with PFC 3D .The sample was formed at zero stress, with weak cementation mimicked by parallel bonds.The linear contact law of the cemented bond acts in parallel with a nonlinear Hertzian contact law.After bond failure, the contact will become purely Hertzian.The figure shows velocities during simulations of (a) isotropic (hydrostatic) stress, (b) uniaxial strain, and (c) a triaxial test.P-wave velocities were measured along the z-and x-axes ("pz" and "px"), and S-wave velocities were measured for propagation along the z-axis with particle motion in the x-direction ("szx"), and for propagation along the x-axis with particle motion in the y-direction ("sxy").
Transmitter Receivers During hydrostatic loading (Fig. 18a), all velocities show a slight increase with stress.The main increase occurs at low stresses, when the Hertzian part of the contact stiffness is most stress sensitive.Since the stress increase is isotropic, the velocity changes are also largely isotropic.The behaviour shown in the figure is strikingly similar to what is seen in many laboratory tests with sandstones (e.g.Holt et al., 1991).The degree of stress sensitivity will depend on the stiffness of the cemented part of the bonds: if that is made higher, then the velocities will become less stress dependent.In the uniaxial strain test (Fig. 18b), the sample actually fails at 6 MPa axial stress, corresponding to a vertical strain of ~ 3 millistrain.The lateral stress at failure is merely 0.4 MPa.At failure, all velocities drop.The velocity reductions are however anisotropic, such that waves propagating perpendicular to the major principal stress ("px" and "sxy") are much more affected than waves propagating parallel to it ("pz" and "szx").This anisotropy is caused by the orientational distribution of bond failures: when bonds break, the stress anisotropy in the material makes it much more favourable for a vertically oriented bond to fail than a horizontal.Bond breakage leads to reduced bond stiffness, and hence reduced velocities.As the confinement is increased, like in the triaxial test shown in Figure 18c (lateral stress = 2.85 MPa), the development of vertical bond failures becomes more suppressed, and the stress-induced velocity anisotropy is reduced.Both these observations show a striking similarity to those made for sandstones in Holt et al. (1991).
The work described above must be considered as preliminary.We have however demonstrated that discrete particle modelling is a promising tool for addressing the stress dependence of wave velocities, an issue that is of central importance when evaluating the feasibility of timelapse seismics in reservoir monitoring studies.As was pointed out by Makse et al. (2001), the discrete particle approach is superior to traditional effective medium theory when it comes to predicting the stress sensitivity of the shear modulus in uncemented grain packs.

CONCLUSIONS
This paper presents some developments and applications with a particle model to study reservoir mechanics.
The behaviours of virgin rock are simulated with a numerical sample formed under in situ stress state.The simulations illustrate the nonlinear stress-strain relationship and stress path of the virgin rock during compaction.Such information is often lost in core tests due to core damage.
Breakable superparticles with irregular shapes are implemented into the particle model and applied to study the failure modes of reservoir rock.The simulations indicate that the failure modes of porous rocks will change with confining stress.With low confining stress, the numerical samples fail with shear band formations.When the confining stress increases, the failure formations tend to become normal to the maximum principal stress.Compaction band is observed when the confining stress is high enough (e.g.18 MPa in a numerical sample with 9.9 MPa UCS as presented in this paper).
A fluid coupling scheme is developed and implemented to 2D and 3D particle models by introducing an associated flow network model.The fluid coupled particle model is applied to study stress-dependent permeability.The results of the simulation fit Kozeny-Carman's Equation.The simulation by 3D model is quantitative, while the 2D is able to do qualitative analysis at this stage.
Wave propagation is simulated with a particle model.It is used to study the acoustic properties of weakly cemented sandstone.The results show qualitative consistence with experiments.
Summing up, the preliminary results of various simulations with a particle model prove the feasibility of applying particle modelling as a numerical laboratory to study particle scale reservoir mechanics.

Figure 4
Figure 4 Development of shear band in particle model at low confining stress.

Figure 5
Figure 5 Development of compaction band in particle model at high confining stress.
Figure 8Flow path.A virtual third particle has to be assumed to calculate the flow path in 2D model.

Figure 9 In
Figure 9 In 3D model, a pore is created between 4 particles.The pore is enclosed in a tetrahedron with 4 particle centres as vertexes.It has 4 pipes connected to its neighbour.Each pipe corresponds to a face of the tetrahedron.(a): a schematic pore; (b): a part of the tetrahedron mesh of a model.
Figure 7 Generation of 2D flow network model.(a): pores firstly generated, enclosed in triangles; (b): cluster pores enclosed in polygons after merging.
Figure 10Setup of the numerical sample under stresses for permeability study.

Figure 11 Pore
Figure 11Pore pressure distribution in a 2D model when the flow is in steady state.Pore pressure is illustrated by disks with their sizes proportional to the values of pore pressures at the locations.(Note: particles are not shown).
Figure 14 (a): stress-strain curves in 2D modelling for permeability alteration and anisotropy study.They are different because the sample dimension in y direction is twice of x direction as in Figure 10; (b): permeability alterations in different loading directions and comparison with estimation by Kozeny-Carman's Equation.
Figure 17 Snap shots of wave fronts.The arrows indicate the velocities of the particles.The figure schematically show the propagating P-wave (a) and S-wave (b).
Figure18(a): acoustic velocities as a function of hydrostatic pressure in a bonded particle model; (b): acoustic velocities as a function of vertical strain during the uniaxial strain test; (c): acoustic velocities as a function of vertical strain during the triaxial loading (constant confining stress 2.85 MPa).