Review of Sand Production PredictionModels

Sand production in oil and gas wells can occur if �uid �ow exceeds a certain threshold governed by factors such as consistency of the reservoir rock, stress state and the type of completion used around the well. e amount of solids can be less than a few grams per cubic meter of reservoir �uid, posing only minor problems, or a substantial amount over a short period of time, resulting in erosion and in some cases �lling and blocking of the wellbore. is paper provides a review of selected approaches and models that have been developed for sanding prediction. Most of these models are based on the continuum assumption, while a few have recently been developed based on discrete element model. Some models are only capable of assessing the conditions that lead to the onset of sanding, while others are capable of making volumetric predictions. Somemodels use analytical formulae, particularly those for estimating the onset of sanding while others use numerical models, particularly in calculating sanding rate. Although major improvements have been achieved in the past decade, sanding tools are still unable to predict the sand mass and the rate of sanding for all �eld problems in a reliable form.


Introduction
A signi�cant proportion of the world oil and gas reserves is contained in weakly consolidated sandstone reservoirs and hence is prone to sand production.Material degradation is a key process leading to sanding.Drilling operations, cyclic effects of shut-in and start-up, operational conditions, reservoir pressure depletion, and strength-weakening effect of water may gradually lead to sandstone degradation around the perforations and boreholes.High pressure gradient due to �uid �ow also facilitates the detachment of sand particles.In addition, �uid �ow is responsible for the transport and production of cohesionless sand particles or detached sand clumps to the wellbore.
Sand production is the cause of many problems in the oil industry and it affects the completion adversely.ese problems include, but are not limited to, plugging the perforations or production liner, wellbore instability, failure of sand control completions [1], collapse of some sections of a horizontal well in unconsolidated formations, environmental effects, additional cost of remedial and cleanup operations, and pipelines and surface facilities erosion, in case the sand gets out of the well.e mechanical prevention of sanding is costly and leads to low productivity/injectivity.erefore, there is always a cost bene�t if sand management and modeling is implemented early before well completions.
Sand production takes place if the material around the cavity is disaggregated and additionally, the operation of the well generates sufficient seepage force to remove the sand grains.It is a complex phenomenon which depends on various parameters such as the stress distribution around the wellbore, the properties of the rock and �uids in the reservoir, and the completion type.erefore, capturing all the factors and mechanisms in the numerical models is difficult and the models have many limitations.
Due to the importance of the sand production prediction in the oil industry, considerable efforts have been made in developing robust numerical methods for sand production prediction.In this paper, the techniques, advances, problems, and the likely future development in numerical models for the prediction of sand production are presented.

Common Techniques Used in Sand Management Decisions
A number of approaches have been developed to predict or help to understand the sand production problem using physical model testing, analytical and empirical relationships, and numerical models.Routine laboratory tests can only predict the onset of sand production [2].More sophisticated physical model could predict volumetric sand production [3].ey are also time-consuming and expensive.In addition, because of the small sizes of the laboratory setup, the results are usually in�uenced by boundary effects.Analytical models are fast and easy to use but they are only suitable to predict the onset of sand production and they have limitations.Most of them are only valid for capturing a single mechanism of sanding and under simpli�ed geometrical and boundary conditions which are not usually the case in complicated �eld-scale problems.Numerical models are by far the most powerful tools for predicting sand production.ey can be combined with analytical correlations to obtain the results more efficiently.Experimental results are also utilized to calibrate or validate a numerical model.Yet, numerical models have their own limitations and extensive efforts have been made to improve them.Modeling of sand production requires coupling of two mechanisms.e �rst mechanism is mechanical instability and degradation around the wellbore and the second one is hydromechanical instability due to �ow-induced pressure gradient on degraded material surrounding the cavity (e.g., perforation and openhole).In general, numerical methods in the mechanical modeling are categorized under continuum and discontinuum approaches.
In the continuum approach, matters are treated as continuous in deriving the governing differential equations.e assumption of continuity implies that the material cannot be separated or broken into smaller pieces.In the case when there is a discontinuity, the magnitudes of deformation along or across the discontinuity are about the same as the rest of the continuum [4].
Discrete element method (DEM) is a useful tool to simulate sand production especially to understand the mechanism of sanding.However, it cannot be used for largescale problems because of large computational time required.e calibration of the model is also difficult and involves several uncertainties as it is not possible to create a model with the exact particle arrangement as the real material.Further, methodologies to directly measure sandstone micro properties have not been developed yet.Presently, the microproperties are obtained in calibrating against the actual sand behavior [5,6].erefore, continuum-based models are more popular especially for �eld-scale problems.However, there can be advanced models which couple continuum and discontinuum models together to take advantage of both methods.ese are known as hybrid models and are discussed later.
A comprehensive sand management may require the use of some or all of the above mentioned techniques.

Numerical Models Based on Continuum Approach.
Developments of continuum models are based on various assumptions, constitutive laws, sanding criteria, and numerical procedures with different levels of complexity to capture the physical behavior of the material.
Table 1 summarizes the majority of continuum-based sanding models.Initially many researchers �rst calculated the onset of sand production or the initiation of mechanical failure until Vardoulakis et al. [7] proposed the basic theory for hydrodynamic erosion of sandstone which is based on �ltration theory without solving the equilibrium equation.Later Papamichos and Stavropoulou [8] combined the evolution of localized deformation with hydrodynamic erosion.is was the beginning for many researches that adopt full strength hardening/soening behavior of sandstone in their models [9][10][11][12][13][14][15][16][17].e results are highly mesh-dependent for strain soening material and hence a regularization method is necessary.Regularization methods include an internal length, which has been related to the grain size, in the formulation.Papanastasiou and Vardoulakis [18] applied Cosserat microstructure method [18,19] for cavity failure around boreholes.Zervos et al. [20] considered Gradient elastoplasticity [21] for thick-walled cylinders.Fracture energy regularization technique [22] was also applied by Nouri et al. [16], Wang et al. [23], and Rahmati et al. [24] in sand production modeling.

Constitutive Models Used in Continuum Approaches.
Rock failure and/or degradation in commonly accepted as prerequisite for sanding.Failure of geomaterials is usually associated with formation of shear bands which are narrow zones of concentrated plastic deformation.is phenomenon, which is known as "deformation localization" or simply "localization, " is one of the key parameters in sanding prediction models.Further details about this concept can be found in Sulem et al. [40], Nouri et al. [16], and Jafarpour et al. [41].
In the simplest form, an elastic brittle failure model has been implemented in sand production models by Nordgren [42], Coates and Denoo [43], Risnes et al. [44], and Edwards et al. [45].Most of these models predict the onset of sand production by considering failure of the sand matrix.Elastic brittle failure rock behavior leads to excessive stress concentrations at the borehole wall and therefore results in overestimation of initial sand production conditions.is model may be used as a quick estimate of sanding onset in relation to production parameters and in situ stresses.e predictions of the elastic models are cumbersome unless they are combined with apparent strength models.
An elastoplastic material model can simulate the material behavior more realistically.However, it requires more computational effort and more input data.Many researchers have implemented elastoplastic models in sand production analysis (e.g., Morita et al. [25]; Antheunis et al. [46]; T 1: Summary of the numerical works on sand production in the literature (continuum approach).[39]).ese works will be discussed in more detail in the next section.Implementation of classical elastoplastic continuum models and the abovementioned failure criteria are inefficient in modeling the localization phenomena, which are of discontinuum nature.erefore, employment of such models in simulation of material degradation leads to inability in recovering size and dependency of the results on numerical mesh design.is problem would be addressed by pursuing two distinct solutions: using discontinuum models that we discuss in the next sections and enriching the material model by an appropriate regularization strategy.Papanastasiou and Zervos [52] used Cosserat continuum and gradient elastoplasticity theories to predict the localization phenomenon.ese models also proved efficient in capturing the size effect oen observed in laboratory tests performed on thick-walled cylinder samples [53].

Geometry
Vardoulakis et al. [54] used bifurcation theory to predict the failure mode and the critical value of the stress at in�nity.ey showed that failure depends on stress path and boundary conditions.Tronvoll et al. [55] performed �nite element modeling using bifurcation theory to solve the axisymmetric problem for the trivial solution and checked the condition for non-axisymmetric bifurcation modes.Van Den Hoek et al. [56] and Papanastasiou and Vardoulakis [57] used bifurcation theory in a Cosserat continuum which takes into account the material microstructure.In a Cosserat continuum individual points possess, in addition to their translational degrees of freedom, independent rotational degrees of freedom, which result in an internal characteristic length in the classical elastoplastic constitutive laws.
e models based on bifurcation theory require special numerical techniques to capture size effects, localization mechanisms, and so forth, which make it computationally demanding and hard to apply in solving �eld problems.
As shown in Table 1, the most basic improvement in sand models is the yield function.Mohr-Coulomb (MC) model is the most common model being used.Vaziri et al. [10] modi�ed the MC model using a bilinear yield function to differentiate sand behavior under low and high con�ning stresses.e model was later used by Nouri et al. [12], Nouri et al. [38], and Vaziri et al. [14].e theory is based on Sulem et al. [40] and is described more thoroughly by Nouri et al. [16] and Jafarpour et al. [41].Haimson and Lee [58] showed that the slit mode of cavity development is related to the formation of compaction bands, which are thin bands of localized compressive deformation in high porosity rocks.erefore, Detournay [15] extended Detournay et al. [11] work by accounting for the compaction mode of failure.e model used Double Yield cap constitutive law to capture compaction bands.e simulation results show that the slit mechanism develops as a combination of volumetric collapse and transport of failed material by seepage forces.It was found out that pore collapse is the responsible mechanism for slit mode of cavity failure associated with sand production.
Although it was considered as one of the main sanding mechanisms, many researchers have avoided the incorporation of compression yielding to simplify their models.
It appears that the optimum constitutive models are those that are based on the critical state theory and use a combined isotropic and kinematic hardening model which allows capturing all kinds of failure (shear, tensile, and compressional).In addition, it would be ideal to capture the effect of hysteresis to simulate fatigue in cyclic start-up and shut-in conditions in the wellbore.

Sanding Criteria Used in Continuum
Approaches.Several mechanisms are recognized as responsible for sand production.ey are mainly based on shear and tensile failure, critical pressure gradient, critical drawdown pressure, critical plastic strain, and erosion criteria.Table 2 summarizes the main sanding criteria used in sand models.
When the effective minimum principal stress is equal to the tensile strength of the formation rock, tensile failure may occur.is mode of failure is responsible for rock degradation.It can occur as a standalone degradation mechanism or in combination with shear failure [22].Tensile mode is also believed to be the responsible mechanism for particle removal aer the degradation during production.In this case, tensile failure relates to seepage forces on the degraded sand particles.
Shear failure may occur when some planes in the vicinity of the wellbore are subjected to higher stress than they can sustain.is is the dominant mechanism in cemented sands and when it is combined with tensile cracks and high compressive stress, it can lead to buckling at the wellbore wall [18].
When effective hydrostatic stresses are increased due to the reservoir pressure depletion, pore collapse may occur which can lead to sand production.Plastic volumetric compression can be captured using a compression yield cap in the failure envelope.is mainly occurs in high porosity sandstones.
Risnes and Bratli [59] proposed a tensile failure criterion for perforation tunnel inner shell collapse.Bratli and Risnes [60] and Weingarten and Perkins [61] proposed sanding criteria in terms of pressure gradient.Morita et al. [25] proposed a sand production model that can be triggered by either shear failure or tensile failure.
Dynamic seepage drag forces lead to internal and surface erosion that result in releasing and transporting sand particles.Internal erosion may be related to micromechanical impacts imposed on solid skeleton by gas bubbles, water droplets, and so forth.Surface erosion may be related to parallel �ow scouring the surface and normal �ow over the surface [28].Numerous authors studied surface erosion criterion in sand production modeling.Vardoulakis et al. [7] proposed a model that was based on mass balance of the produced solids and radial �owing �ow conditions, the constitutive law for particle erosion, and Darcy's law but neglecting skeleton deformation.Based on sand production experiments, Tronvoll et al. [62] showed that in addition to the radial �ow, axial �ow parallel to the perforation is can be considered as an elaboration of tensile criterion considering the effect of friction coefficient e element is removed suddenly without allowing the porosity to grow gradually channels is important in sand production and it may result in surface erosion of the perforation channels.Consequently, Vardoulakis et al. [63] extended Vardoulakis et al. [7] work to account for axial �ow conditions.In the governing equations, they included Brinkman's extension of Darcy's law, which accounts for a smooth transition between channel �ow and Darcian �ow.e results show that erosion progresses in time at the front of high transport concentration.
First Stavropoulou et al. [64] and later Papamichos et al. [9] advanced a purely hydromechanical model proposed by Vardoulakis et al. [7] by coupling the poromechanical behavior of the solid �uid system with the erosion behavior of the solids due to �uid �ow.Papamichos et al. [65] extended their own work by using a porosity diffusion type law that results in a sand rate that decreases over time when the process of erosion zone enlargement takes place.e model is based on nonlinear elastoplasticity, nonlinear stress dependent elasticity, friction hardening and cohesion soening, and single-phase �ow fully coupled with geomechanics.
Wang et al. [35] also performed a fully coupled singlephase �ow analysis based on erosion mechanics using the FEM.ey applied the model in 2D plain strain geometry for both open hole and perforated casing completion.
Based on laboratory experiments, Haimson [66] and Papamichos [67] observed the slit cavity evolution pattern during sand production.Subsequently Detournay et al. [11] proposed a model to predict the formation of slit channels based on a critical �ow rate.ey modi�ed the erosion law used by Vardoulakis et al. [7] with the addition of a critical �ow rate which is a function of the grain size.ey also assumed that sand continuously is produced until a critical porosity is reached beyond which the material suddenly collapses.is process can be responsible for the periodic sand bursts observed in experiments.e model was applied to 2D plain strain geometry for a long wellbore or perforation using a �nite difference soware.e model can predict different erosion features such as surface spalling and small burst events.
Kim et al. [17] calculated sanding conditions and utilized a sanding criterion based on the force balance on each element and achieved a good match with experimental data reported by Nouri et al. [51].In their criterion, the forces leading to sand production are hydrodynamic forces generated by the pore pressure gradient between element faces in the �ow direction.e resistance forces are the forces that retain the elements in place and are generated by vertical and tangential stresses and the friction coefficient.e advantage of using this method is that no calibration parameter for sanding (such as sand production coefficient) is necessary.e friction coefficient is an empirical parameter which depends on the grain size and mineralogy.
However, the abovementioned coupled hydromechanical erosion models have been criticized due to the following con�icting assumptions.Material mass balance equations are based on rigid porous media while equilibrium equations are based on deformable porous media.erefore in order to establish a proper coupled mechanical erosion model in a consistent manner, the porosity changes can be split into two parts: one related to volume changes as a result of erosion, and the other one due to deformation in the matrix subjected to stress changes [48].
To sum up, a realistic sanding model should ideally account for all failure mechanisms (shear, tensile, and compressional) and must also consider the effect of �uid �ow.erefore, a suitable sand erosion model consists of a combination of erosion criterion, tensile criterion, and compressional criterion.Considering the physics of sand production, erosion seems more suitable for weak rock where full decementation and degradation to small particles is more likely [68].On the other hand strong rocks are more prone to localized failure that result in larger chunks of sandstone which are not easily eroded away.Lastly, compressional failure is more dominant in highly porous weak materials where void spaces collapse easily under high loading.

Phases Involved in Continuum
Approaches. e models can be categorized into two groups based on the phases involved.In the �rst group, mass balance equation is solved for only �uid and solid phases while the second group recognizes �uidized solid as a phase and solves for solid concentration in the �uid.Fluidized solids are the particles in suspension that move with the �uid.Any other loose particle which is trapped inside the void space is seen as part of the solid phase.However, these models use a constant viscosity for the slurry.It is notable over time that researches tend to use the simpler approach (solid and �uid phases) and couple the equations with an erosion model.is is mainly because good agreement between the model and �eld observations were obtained when combined with a suitable sanding criterion.
Multiphase �uid �ow may also affect both the onset of sand production and sand rate.Tronvoll et al. [69] and Vaziri et al. [10] observed water cut effects on the onset of sand production.Water in�ow changes the relative permeability and capillary pressure.It also can dissolve cement bonds and weaken the strength of the porous media.
Wang et al. [48] presented an integrated modular approach to predict volumetric sand production and cavity growth under two-phase �ow (oil and water) and 3D geometry.In this work, the effect of water on rock strength reduction is re�ected in material properties such as cohesion.e results show that water contact has a signi�cant impact on the sand rate.
�as �ow also may hasten the instability process in sand production.When gas comes out the oil phase due to pressure drop and �ows towards the wellbore at high velocities, it applies additional drag forces on sand particles and increases sand production.Wan and Wang [49] studied the gas effects in 1D model using the �nite element method.ey assumed that eroded mass in erosion law is proportional to the total �uid �ux, which is referred to the oil and gas �ux.e effect of multiphase �uid �ow was also considered by Wang and Xue [32] for oil-water phases and later by Wang et al. [35] for oilwater-gas phases.
e only available numerical works in the literature that incorporate the water contact effect on sanding were performed by adjusting the cohesion or lowering the rock strength [17,39].

2.1.4.
Treating the Sanded Elements.Different strategies have been used for dealing with those elements that satisfy sanding criteria.e �rst one is to remove such elements from the model assuming cavities and wormholes grow as a result of sanding.is seems to be a suitable approach in stronger rocks in which stable cavities can form.e other approach is to keep the element in place but alter the material properties to residual values.We believe that stable cavities cannot develop in weak sandstone.Rather, the space is occupied by in�ll material, or cohesionless sand particles.In this approach, sandstone properties are altered to those of degraded cohesionless sand.e property change should also be applied in erosion models as a function of increasing porosity in the eroding elements.
It is evident that the moduli, tension cut off and permeability vary with the production of sand and increase of the porosity.However, the correct method to apply these changes requires experimental data.Most researchers use arbitrary choice of permeability change with porosity or with volumetric strain or even with mean stress.Wang and Xue [32] used two permeability correlations and found that the permeability relationship plays an important role.ere is also disagreement on whether permeability should increase or decrease.It is reported that for high permeable sand formations, the permeability of the disaggregated sand is much less than that of the intact sand.is is mainly attributed to the sand deposition and plugging of the pore space, which is not the case for less permeable sand [26].
e moduli of the elements will also vary with porosity to the values of loose cohesionless sand (about 6.9 MPa for bulk modulus and 4.14 MPa for shear modulus).ese numbers are the lowest values reported for loose sand [70].

Model Design.
Openhole completion is oen treated with axisymmetrical models.Strictly speaking, this is correct only when horizontal stresses are equal.However, most of the times, if not always, principal horizontal stresses,  ℎ and   , are not equal.In such cases a plane strain model can be a suitable choice for 2D analysis.Plane strain may not always be an appropriate assumption as vertical deformation may not necessarily be negligible.Papanastasiou and Zervos [52] suggested that generalized plane strain could be an appropriate assumption in the modeling of vertical and inclined wellbores.
As sand production model is commonly used in cased and perforated (C&P) completions, it is important to consider the sand behavior under such geometrical and boundary conditions.For instance the frequency of shots, the length of perforations, and their orientations may lead to a more intense commingling of the failed zones and �nally higher sanding rate.One solution is using 3D simulation but it is computationally demanding because very �ne mesh is required around perforations.In addition, creating the hollow geometry in weak rocks may not be reasonable as in reality the perforation can collapse upon creation and be �lled with in�ll degraded sand.
Some perforation simulations have used axisymmetric models in which the perforation is assumed to be ring-shaped rather than conical or cylindrical [39].Such an assumption in�uences the pressure gradient around the perforations and also impacts the mechanical response.ese models are also unable to capture the perforation direction effect, which can play a signi�cant role in perforation stability and failure as demonstrated by Papanastasiou and Zervos [71].ey performed a 3D numerical simulation to study the effect of orientation on stability and failure of perforation.Based on results of this work, it is recommended to avoid perforating the wellbore parallel to the minimum horizontal stress direction as perforations in this direction suffer more compressive stress and hence more chance of failure and sand production.ese models require calibration against �eld and�or physical model testing before application to realworld sanding problems.

Other Factors for Continuum-Based Numerical Models.
Sand production is a moving-boundary problem.As sand is produced, a sanded zone is formed around the perforations.Adaptive meshing can be very useful in processes where the geometry or the boundary is changing.Yet, there are only two applications of adaptive meshing in the literature [23,37].
e current models are unable of predicting the sand bridging and �nes retention in the rock.Sand particles can aggregate at the perforation cavity and form a stable entity called sand bridge and act as a �lter which may reduce further sand production as long as the �ow rate remains constant.e theory for modeling sand deposition was proposed by Vardoulakis et al. [7] but it has been applied only in one model [30] using a similar but not exactly the same approach.Yi [30] considered a part of the degraded sand as to be deposited in the porous media.e difficult part about modeling sand deposition is the calibration of the critical porosity or the critical sand concentration aer which sand deposition initiates.
An important issue about the current sand models is that almost all of them are applied in modeling production wells.A few researchers [14,30] performed numerical studies to predict sanding in injector wells.Sanding mechanisms in injectors have not been investigated thoroughly and may be quite different such as the effects of water hammer (WH) waves.e observation in injection wells is oen described as the cases with high sand production within a short duration.Few papers [72][73][74][75] tried to explain the sanding problems in injectors.It is stated that sand liquefaction due to WH pressure pulses is the most likely mechanism for massive sand production.Liquefaction is de�ned as the process by which saturated sand loses shear strength and stiffness in response to dynamic loading [76].WH is a general term describing generation, propagation, and damping of pressure waves in pipes.It occurs due to sudden velocity changes such as quick shutting of the well [72].Nevertheless, no work has been published which investigates liquefaction around the wellbore and the conditions leading to liquefaction.As a result, sand particles can �ow easily like a liquid.As no investigation has been performed on liquefaction in sand production, it is difficult to con�rm its role in massive sand production.

Numerical Models Based on Discontinuum Approach.
Sand production is a continuous and dynamic process that occurs at the microscopic scale and the rock becomes a discontinuum in nature.As mentioned before, conventional continuum approaches cannot capture local discontinuous phenomena.erefore, discontinuum approach is promising to simulate phenomena such as detachment of individual particles from the rock matrix.
Cundall [83] �rst introduced the Discrete Element Method (DEM).e method can be used to simulate the disintegration of granular media subjected to loading.Each particle of the granular media is considered as an individual entity with a geometric representation of its surface topology and a description of its physical state.Particle bonds are modeled with a spring-dashpot in the normal direction and a spring-dashpot-frictional slider in the tangential direction.In the DEM, the interaction of the particles is treated as a dynamic process and a state of equilibrium is reached whenever the internal forces are equal to the external forces.e contact forces and displacements of a stressed assembly of particles are found by tracing the movements of the individual particles [84].
Table 3 summarizes some of the discontinuum-based sanding models.At �rst, ��Connor et al. [77] introduced the application of DEM to model the mechanics of sand production during oil recovery.Using laser scanning and sieve testing they developed techniques to consistently represent particles with irregular geometries.ey used particle bonding scheme to mimic the cement and cohesion due to capillary forces.e bond also incorporated spring stiffness and a nominal tensile breaking strength assuming a bond dimension proportional to the size of the particles it connects.ey incorporated the �uid �ow calculations by combining the continuity equation and Darcy�s law using the �nite element method.Darcy �ow is formulated with a measure for effective permeability in the solid medium based on the porosity and average diameter of the solid particles.eir 2D model provides an understanding of the fundamental physics involved in sand production and the relative importance of various rock and �uid properties.
Jensen and Preece [78] explored the coupling of 2D DEM and �nite element implementation of the 2D continuity equation for Darcy �ow to assess the sanding potential.e basic particle shape used by the model was an nsided polygon and only the tensile mode for bond failure was considered.ey concluded that lower strength of the cohesive bonds increased the number of particles breaking free from the solid matrix.
Li et al. [79] used commercial DEM code PFC2D to simulate hollow cylinder tests with �uid �ow to study sanding.PFC2D simulates an assembly of circular disks with the bonds inserted between them.e disks are assumed to be rigid but they can overlap.e bonds have normal and shear stiffness and strength.In the standard PFC2D code, a bond fails when the tensile or shear stress in the bond exceeds its strength.Bond breakage may be interpreted as microfailure in the real rock.e growth of such microfailures eventually leads to macroscopic failure of the rock.Li and Holt [80] showed that DEM model may not result in realistic macroscopic friction coefficients if only circular or spherical grain shapes are used.
Li et al. [79] improved the prediction of macroscopic friction coefficient by setting the bond strength so high that no bonds in the model would fail due to the stress in the bond.Instead, all bonds associated with a given disk break when the stresses inside the disk satisfy a failure criterion.ey used a failure criterion composed of tensile failure, shear failure, and compressive failure.A simple approach was T 3: Summary of the numerical works on sand production in the literature (Discontinuum approach).

Model
Geometry used to calculate and couple the �uid �ow.ey found three typical failure patterns in the simulations similar to those observed in laboratory experiments.e slit-like breakout failure pattern was observed when the material is prone to localized compressive failure due to grain crushing.For those cases where the material was weak and tensile strength was low, uniform failure around the borehole was observed along with a rather uniform hole enlargement.In those cases with relatively competent rock properties, which were unlikely to fail in localized compaction, the failure pattern was observed to be in the form of dog-eared breakouts.
Several researchers (e.g., [80,85]) have modeled �uid �ow in 2D DEM codes by introducing �uid �ow networks and simulating �ow along the �ow paths connecting the voids which are referred to as pipes.e �uid velocities that �ow through the pipes and the pore pressures were computed based on Darcy's theory.e forces arising from the pore �uid were then calculated and applied to the particles in the DEM model.Li and Holt [80] implemented this type of �uid-solid coupling system in the PFC2D codes.�hile the geometrical limitations of using a 2D model to investigate breakout geometry are obvious, the �uid �ow simulation through �ow networks is computationally expensive.e entire coupling techniques described above used Darcy's law to calculate �uid �ux or pressure.Darcy's law has been derived from the Navier-Stokes equations via homogenization that is only valid for slow and viscous �ow.e conditions may not be met around the wellbore where breakout forms [81].Chan and Tipthavonnukul [86] proposed a method to couple continuum and discontinuum �ow to simulate granular particles movement in a �owing �uid.e �uid �ow is modeled using the 2D Navier-Stokes equations solved by a �nite volume method.e coupling is achieved by detecting the presence of the solid in the �ow domain and altering the �ow resistance accordingly.e method is computationally expensive and therefore it is not applicable in large-scale problems.
In 2D DEM models, only two force components and one moment component are considered.However three force components and three moment components exist in a 3D DEM model.Since sand production is a three-dimensional problem with complex geometry, where cement bonds and arching of particles are important, there is a need to model this problem using 3D DEM.e 2D DEM models overestimate the effect of �uid �ow on the integrity of the assembly of particles as the resistance to particle dislodgment due to contact forces and bonds normal to the �uid �ow is neglected in 2D models.In addition, 2D models cannot represent threedimensional pore �ow networks.
Cheung [81] used a coupled �uid-solid 3D DEM model for a perforation test simulation to study the sand production problem.ey used the Navier-stokes equations assuming radial �ow.Although this scheme is computationally inexpensive, it has two major problems.First, by assuming radial �ow, it is not suitable for investigating the impact of the �ow at the tip of the perforation where the �uid �ow is in all directions.Second, the �uid �ow scheme does not consider the presence of the cement between particles.e magnitude of the radial pore �uid velocity in each �uid cell is calculated, considering only the presence of the particles.e presence of cement can highly affect the rock conductivity and therefore the �uid velocity.
Later, Zhou et al. [82] employed DEM with computational �uid dynamics (CFD) and showed that the main features of sand erosion can be captured by the CFD-DEM approach.
e main advantage of the DEM models is that it captures the motion and interaction of individual sand grains and its failure micro mechanism in a dynamic process.It enables the model to predict many real behaviors such as continuously nonlinear stress strain response, behavior that changes in character according to stress state, memory of previous stress or strain excursion in both magnitude and direction, dilatancy that depends on history, mean stress, and initial states, hysteresis at loading/unloading among others.
To the best of our knowledge, there is no existing continuum constitutive model that reproduces all of these behaviors.However, since the DEM involves many individual particles and interactions between them, it is computationally expensive and therefore it is not applicable to large-scale problems.
Another disadvantage of the DEM model is the lack of a systematic method for an objective determination of micro material parameters.As opposed to the continuum-based models for which the strength and elastic properties can be determined directly from laboratory testing, the micro properties cannot be determined by direct measurements of the macro responses on the laboratory specimens.It could be found by means of a calibration process in which a particular assembly of particles with a set of micro parameters is used to simulate a set of material tests and the micro parameters then are evaluated to reproduce the macro responses measured in such tests [84].Several researchers (e.g., [5,6,84,87]) have proposed calibration procedures relating the micro parameters to macro properties of the material.However, the calibration procedure is challenging for the 3D DEM models.ere may be several variations in the parameters and it is difficult to conclude which set of parameters is most appropriate for the material.

Hybrid Approaches.
Considering the advantages and disadvantages of the continuum-and discontinuum-based approaches, a hybrid model that combines them can be practical and efficient in sand production modeling.Continuumbased approaches can be used in far �eld where the deformation is small hence continuum assumption is still valid and computationally efficient.On the other hand, a discontinuum-based approach can be used to describe large deformation or discontinuity near the wellbore or the perforations.In this manner, accurate and descriptive simulation of �eld-scale problems becomes possible with nowadays computational power available.
Some researchers have used this approach to analyze geomechanical problems.For example, El Shamy and Elmekati [88] and Elmekati and Shamy [89] combined a FEM code with a DEM code to analyze soil-structure interaction problems.Also, Azevedo and Lemos [90] used the same approach to study fracture growth in tensioned columns.In a similar work, Zeghal and El Shamy [91] coupled a continuum �uid model with discontinuum particle model to analyze the dynamic liquefaction of granular soils.
To the best of our knowledge, the hybrid scheme has not been used in sand production modeling.

Conclusions
Despite the numerous efforts in sand production and modeling, there are still some fundamental de�ciencies which require to be addressed.Considering the works reported in the literature, there is still signi�cant room for improvement in sanding models.Some are listed below.
(i) Critical-state-based constitutive models are expected to provide a more realistic representation of sand production.Such models can also predict compressive yielding which is the main mechanism for the creation of slit mode of degradation that may precede sanding.Yet, researchers have tended to avoid them due, perhaps, to several calibration parameters.
(ii) Sanding criteria which are based on erosion mechanics are most popular among researchers for describing and simulating the sanding phenomenon.ese criteria usually involve a constitutive law for mass generation that requires calibration against laboratory tests.e erosion criteria are then coupled with hydromechanical behavior of the �uid-rock system for the purpose of numerical simulation.Erosion sanding criteria are more applicable for weak rocks where degradation of the rock into small particles is more probable.For strong rocks, other sanding criteria such as those that are based on tensile or shear failure or their combination with the erosion criteria could be more appropriate.
(iii) e most appropriate constitutive law may be a combined isotropic and kinematic hardening model, which allows capturing of the principal mechanisms of failure (shear, tensile, and compressional) as well as capturing the hysteresis effect in multiple shut-in and start-up scenarios.
(iv) Simulation of perforated wellbores requires a very �ne mesh around perforations which makes the analysis computationally demanding especially is 3D modeling.2D models of perforation using plane strain and axisymmetric assumptions do not fully capture the accurate mechanical and �uid �ow responses.Yet, the majority of present models are simulated in 2D.A possible solution may be developing special numerical elements that best represent the perforation geometry and response.
(v) In capturing sanding in weak rocks, it is more realistic to keep the elements that satisfy sanding criteria in place and assign residual properties than to remove them.Element removal makes more sense when treating stronger rocks.
(vi) e assessment of the properties of in�ll materials requires more investigations and experimental data.e choice of these properties plays a signi�cant role in sanding prediction.More accurate correlations for the changes of rock and �ow properties with sand production and increase of porosity are essential.
(vii) Methods to capture sand arching have not been explicitly developed as such require complex interaction between the geometry of the opening in the completion and the disaggregated rock mass characteristics under prevailing stress state.Sand particles aggregate at the perforation cavity and form a stable entity called sand bridge and acts as a �lter that prevents further sand particles to produce as long as the �ow rate remains constant.
(viii) Sand liquefaction around injection wells has not been investigated fully yet.Since there is no recorded measurement of how sanding occurs in an injector, it is not yet clear whether it is occurs suddenly or gradually over many cycles, and if it is due to waterhammer or cross and back �ow, or whether the produced material is sand or shale/clay.
(ix) e calibration procedure in the DEM model for determining micro material parameters needs further research.
(x) In order to have more accurate analysis in the DEM model, the �uid �ow scheme needs to be modi�ed.For instance, in current models, the permeability is usually related to porosity changes due to particles removal.However, cement debonding and wash out could also affect the pore-network and therefore the permeability.
(xi) Hybrid model combining DEM for the rocks around the wellbore and continuum approach for far-�eld rocks is expected to provide a more realistic and yet practical representation of a number of critical factors governing sanding.�on��ct of �nterests e authors of this paper are not involved with those soware companies whose products have been mentioned in this paper that may give rise to con�ict of interests.