# Review of Sand Production Prediction Models

Authors: Hossein Rahmati; Mahshid Jafarpour; Saman Azadbakht; Alireza Nouri; Hans Vaziri; Dave Chan; Yuxing Xiao
Journal: Journal of Petroleum Engineering
 (2013)

---

## Abstract

Sand production in oil and gas wells can occur if fluid flow 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. The amount of solids can be less than a few grams per cubic meter of reservoir fluid, posing only minor problems, or a substantial amount over a short period of time, resulting in erosion and in some cases filling and blocking of the wellbore. This 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. Some models 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 field problems in a reliable form.


---

## Body

## 1. Introduction

A significant 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 fluid flow also facilitates the detachment of sand particles. In addition, fluid flow 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. These 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 clean-up operations, and pipelines and surface facilities erosion, in case the sand gets out of the well. The mechanical prevention of sanding is costly and leads to low productivity/injectivity. Therefore, there is always a cost benefit 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 fluids in the reservoir, and the completion type. Therefore, 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.


## 2. 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]. They are also time-consuming and expensive. In addition, because of the small sizes of the laboratory setup, the results are usually influenced 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 simplified geometrical and boundary conditions which are not usually the case in complicated field-scale problems. Numerical models are by far the most powerful tools for predicting sand production. They 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. The first mechanism is mechanical instability and degradation around the wellbore and the second one is hydromechanical instability due to flow-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. The 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 large-scale problems because of large computational time required. The 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]. Therefore, continuum-based models are more popular especially for field-scale problems. However, there can be advanced models which couple continuum and discontinuum models together to take advantage of both methods. These 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.


### 2.1. 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.Table1 summarizes the majority of continuum-based sanding models. Initially many researchers first 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 filtration theory without solving the equilibrium equation. Later Papamichos and Stavropoulou [8] combined the evolution of localized deformation with hydrodynamic erosion. This was the beginning for many researches that adopt full strength hardening/softening behavior of sandstone in their models [9–17]. The results are highly mesh-dependent for strain softening 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 micro-structure 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.Table 1
Summary of the numerical works on sand production in the literature (continuum approach).
Model
Geometry and solution method
Yield
Hardening/softening
Coupling
Phases
Sanding criteria
Permeability alteration in the sanded zone
Other features
Morita et al. [25], Morita et al. [26], Burton et al. [27]
3D; finite element (FE)(SAND3D software)
Kinematic model with a cap
yes (flow friction)
Iteratively  coupled
(1) Fluid(2) Solid
Maximum plastic strain limit
No change
(1) Only the onset of sanding(2) Burton applied it for gas reservoirs
(Vardoulakis et al. [7])
1D; finite difference (FD)
N/A
N/A
Fluid flow and erosion are coupled
(1) Fluid(2) Fluidized    solid
Erosion
k
/
k
0
=
φ
3
/
(
1
-
φ
2
)  (Carman-Kozeny)
(1) Only hydromechanical effects; equilibrium eqn. is not solved(2) Sand deposition is neglected in modeling
Skjaerstein et al. [28]
1D; FD
N/A
N/A
Fluid flow and erosion are coupled.
(1) Fluid(2) Fluidized solid
Erosion
1
/
k
=
B
φ
+
A  (found experimentally)
Forcheimer’s law was used instead of Darcy’s law to account for turbulence
Papamichos and Malmanger  [29]
2D Axial symmetry and 3D; FE;  Newoton-Paphson (NR) iterations
MC
Yes
Fully coupled
(1) Fluid(2) Solid(3) Fluidized solid
Erosion
k
/
k
0
=
φ
3
/
ρ
-
(
1
-
φ
2
)  ρ-=ρf+c(ρs-ρf)
Tension cut-off: function of both plastic strain and porosity; by the factor(1-φ)/(1-φi)
Papamichos et al. [9]
2D Axial symmetry;  FEM; NR iterations
MC
Yes
Fully coupled
(1) fluid(2) solid
Erosion
Carman-Kozeny
Tension cut-off and Young mod changed by the factor(1-φ)/(1-φi)
Yi [30]
2D axial symmetry;  FD
MC
No
(1) Fluid(2) Solid(3) Fluidized sand
Erosion
k
/
k
0
=
γ
(
φ
3
(
1
-
φ
0
2
)
/
φ
0
3
(
1
-
φ
2
)
)
×
(
1
/
(
1
+
β
m
d
/
ρ
s
)
)
Sand deposition in porous media is considered
Vaziri et al. [31]
2D axial symmetry; FE; fully-implicit
Modified MC with tensile failure
Yes
Fully coupled
(1) Fluid(2) Solid
Tensile failure
k
/
k
0
=
exp
[
-
2.88
×
10
-
3
(
σ
i
′
-
σ
0
′
)
]
Zero stiffness, compressibility and highk for the liquefied tensile zone
Wang and Xue [32]
FE; Crank-Nicholson for time integration
MC
No
Fully coupled
(1) Oil(2) Water(3) Solid(4) Fluidized sand
Erosion
Kozeny-Poiseuille law and Carman-Kozeny
Chin and Ramos [33]
2D and 3D; FE; explicit; NR iteration
Drucker-Prager
No
Coupled
(1) Fluid(2) Solid
Shear dilation
Power law with porosity (exponent = 5.6)
Porosity is changed as a function of plastic volumetric strain
Nouri et al. [34]
2D plane strain;  FD
Bilinear MC
Yes
Fully coupled
(1) Fluid(2) Solid
Tensile failure or shear-failed element falls in tension
0
Capillary is considered as a residual cohesion
Wang et al. [35]
2D; FE
Drucker-Prager
No
Fully coupled
(1) Solid(2) Fluidized solid  (3) Oil(4) Water(5) Gas
k
=
k
0
exp
[
A
(
(
φ
-
φ
0
)
/
(
φ
max
-
φ
0
)
)
]
Cohesion and friction drop linearly with porosity
Servant et al. [36]
2D; FE
MC
No
Iteratively  coupled
(1) Fluid(2) Solid(3) Fluidized solid
Yielding
Not mentioned
Failed material is treated as a Poiseuille fluid.Constant viscosity for the slurry
Detournay et al. [11]
2D; FD
MC
Yes
Iteratively  coupled
(1) Fluid(2) Solid
Erosion
k
/
k
0
=
φ
3
(
1
-
φ
0
2
)
/
φ
0
3
(
1
-
φ
2
)
Bulk mod. Change byφ0/φ
Nouri et al. [12]
2D axial symmetry; FD
Bilinear MC
Yes
Iteratively  coupled
(1) Fluid(2) Solid
Tensile failure or a shear-failed element falls in tension
k
=
C
φ
3
/
γ
f
(
1
-
φ
2
)
Nouri et al. [37]
FE; NR iterations
MC
Yes
Fully coupled
(1) Fluid(2) Solid
Tension
0
Adaptive mesh is used
Nouri et al. [38], Vaziri et al. [14]
2D plane strain and axial symmetry; FD
Bilinear MC
Yes
Iteratively  coupled
(1) Fluid(2) Solid
Complete degradation and tensile mean effective stress
High permeability is assigned to infill materials (elements that satisfy sanding criteria)
WH pulsing is included in the model, and stiffness changes with sanding (Vaziri et al. [14])
Detournay [15]
2D plane strain;  FD
Double-yield
Yes
Iteratively  coupled
(1) Fluid(2) Solid
Same as (Detournay et al. [11])
Not mentioned
Only cap yielding is considered
Nouri et al. [16]
All the features are the same as (Nouri et al. [38]) but fracture energy regularization is applied to remove mesh dependency and the calibration of hardening/softening behavior is discussed thoroughly
Kim et al. [17]
3D;  FD
MC
Yes
Not mentioned
(1) Fluid(2) Solid
Calculated from force balance on the element
N/A (element removal)
No calibration parameter to match the experiment is used and yet a good match is observed
Wang et al. [23]
2D;  FE(ELFEN software)
Soft rock model
Yes
Fully coupled
(1) Fluid(2) Solid
Displacement  (the materials passing the wellface are removed)
K changes as a function of dilation or volumetric strain: k=exp(Aεv)
(1) Adaptive mesh(2) Arbitrary Largrangian Eulerian formulation(3) Fracture Energy Regularization
Azadbakht et al.[39]
2D, FE for fluid and FD for solid
Bilinear MC
Yes
Sequentially coupled
(1) Fluid(2) Solid
Erosion
K changes linearly with porosity
Water effect is considered by reducing the cohesionBulk moduli change linearly with porosityAfter critical porosity, residual properties are assigned to the element as infill material


#### 2.1.1. 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. This 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. This model may be used as a quick estimate of sanding onset in relation to production parameters and in situ stresses. The 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]; Peden et al. [47]; Papamichos and Vardoulakis [9]; Detournay et al. [11]; Wan and Wang [48]; Servant et al. [36]; Wang et al. [35]; Wan et al. [49]; Detournay [15]; Wan and Wang [50]; Vaziri et al. [14]; Nouri et al. [13]; Nouri et al. [51]; Nouri et al. [34]; Rahmati et al. [24]; Azadbakht et al. [39]). These 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. Therefore, employment of such models in simulation of material degradation leads to inability in recovering size and dependency of the results on numerical mesh design. This 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. These models also proved efficient in capturing the size effect often observed in laboratory tests performed on thick-walled cylinder samples [53].Vardoulakis et al. [54] used bifurcation theory to predict the failure mode and the critical value of the stress at infinity. They showed that failure depends on stress path and boundary conditions. Tronvoll et al. [55] performed finite 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.The 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 field problems.As shown in Table1, 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] modified the MC model using a bilinear yield function to differentiate sand behavior under low and high confining stresses. The model was later used by Nouri et al. [12], Nouri et al. [38], and Vaziri et al. [14]. The 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. Therefore, Detournay [15] extended Detournay  et al. [11] work by accounting for the compaction mode of failure. The model used Double Yield cap constitutive law to capture compaction bands. The 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.


#### 2.1.2. Sanding Criteria Used in Continuum Approaches

Several mechanisms are recognized as responsible for sand production. They are mainly based on shear and tensile failure, critical pressure gradient, critical drawdown pressure, critical plastic strain, and erosion criteria. Table2 summarizes the main sanding criteria used in sand models.Table 2
Sand criteria used in continuum-based sanding models.
Sanding criterion
Formulation
Further developments
Assumptions
Vardoulakis et al. [7]
d
φ
/
d
t
=
m
˙
/
ρ
s
=
λ
(
1
-
φ
)
c
∥
q
i
∥  And with sand deposition:= λ(1-φ)(c-c2/ccr)∥qi/q0∥n
Sand deposition is also considered in some models [28]
It is assumed that the sand in place is fully degraded from the beginning and the production is only due to the hydrodynamic forces. Equilibrium eqn. for solid phase is often ignoredTo initiate the process a very small solid concentration is given as a boundary condition. The results are insensitive to this value as long as it is small
Skjaerstein et al. [28]
d
φ
/
d
t
=
m
˙
/
ρ
s
=
λ
c
(
φ
-
φ
2
/
φ
cr
)
q
x
Papamichos and Malmanger [29]
m
˙
/
ρ
s
=
λ
(
1
-
φ
)
c
q
i
q
i  λ=λ(εp) if εp≥εppeak
The critical flow rate is ignoredSoftening is not simulated
Yi [30]
∂
φ
/
∂
t
=
η
λ
(
1
-
φ
-
(
m
d
/
ρ
s
)
)
v
2
-
λ
d
v
c
s  md˙/ρs=λdvcs+(1-η)(m˙/ρs)
Papamichos et al. [9]
ε
p
>
ε
p
cr  m˙/ρs=λ(1-φ)∥qi∥if φ<φcr, else pore collapse or sudden element removal.λ=λ(εp)
Later models definedεpcr as a function of stresses (Wang and Xue [32]). Critical plastic strain εpcr  was later linked to the state of zero cohesion (Azadbakht et al. [39])
Porosity is assumed to increase until it reaches unity. Usually, it is assumed that eroded elements collapse once their porosities reach a critical level less than unity (say 0.5)
(1) Erosioncriterion
Chin and Ramos [33]
d
φ
/
d
t
=
(
1
-
φ
)
∇
·
v
s
Wang [32, 35]
d
φ
/
d
t
=
m
˙
/
ρ
s
=
λ
(
1
-
φ
)
c
q
i
q
i
Wang et al. [35] assumed λ=λ(εp)
The resisting effect of cohesion was ignored resulting in conservative sanding assessmentSoftening is not considered
Servant et al. [36]
After failure, new material properties are assigned to the failed material (such as zero cohesion). The infill is assumed to behave as a Poiseuille fluid with a viscosity depending on the sand concentration
Detournay et al. [11], Detournay [15]
q
>
q
cr  qcr=qcr(k,φ,Rp,∅,C)  if φ<φcr,∂φ/∂t=λ(1-φ)(q-qcr)(SL/VL)  else element removal
Azadbakht et al. [39]
ε
p
≥
ε
p
cr
and
q
>
q
cr  if φ<φcr  m˙/ρs=λ(1-φ)(qi-qcr)λ=λ(εp)qcr=qcr(k,φ,Rp,∅)else pore collapse or element removal
Some models consider sand production from elements whose cohesion is already reduced but not equal to zero yet (Detournay et al. [11]): qcr=qcr(k,φ,Rp,C,∅)
(2) Tensilecriterion
Vaziri et al. [10], Nouri et al. [12],  Nouri et al. [13], Vaziri et al. [14]
Complete degradation (c=0) and tensile failure
Tensile failure has been replaced by either tensile effective stress or tensile mean effective stress, assuming cohesionless sand has no tensile strengthA more comprehensive one is: shear-failed element with zero cohesion falls in tension or sand fails in tension (Nouri et al. [12])
The element that satisfies sanding criteria is removed from the mesh (assuming cavity can grow adjacent to the wellbore)
(3) Force balanceon elements
Kim et al. [17]
A
r
p
∣
r
=
L
-
A
r
∣
r
=
0
p
w
>
A
r
σ
r
∣
r
=
L
r
=
0
-
A
θ
σ
θ
r
∣
θ
=
0
θ
=
L
-
A
z
σ
z
r
∣
z
=
0
z
=
L
+
μ
f
(
A
θ
σ
θ
∣
θ
=
0
+
A
θ
σ
θ
∣
θ
=
L
+
A
z
σ
z
∣
z
=
0
+
A
z
σ
z
∣
z
=
L
)
For illustration purposes the boundary is shown from 0 to L in all directions.This can be considered as an elaboration of tensile criterion considering the effect of friction coefficient
The element is removed suddenly without allowing the porosity to grow graduallyWhen the effective minimum principal stress is equal to the tensile strength of the formation rock, tensile failure may occur. This 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 after 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. This 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. This 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 flow scouring the surface and normal flow 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 flowing flow 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 flow, axial flow parallel to the perforation 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 flow conditions. In the governing equations, they included Brinkman’s extension of Darcy’s law, which accounts for a smooth transition between channel flow and Darcian flow. The 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 fluid system with the erosion behavior of the solids due to fluid flow. 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. The model is based on nonlinear elastoplasticity, nonlinear stress dependent elasticity, friction hardening and cohesion softening, and single-phase flow fully coupled with geomechanics.Wang et al. [35] also performed a fully coupled single-phase flow analysis based on erosion mechanics using the FEM. They 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 flow rate. They modified the erosion law used by Vardoulakis et al. [7] with the addition of a critical flow rate which is a function of the grain size. They also assumed that sand continuously is produced until a critical porosity is reached beyond which the material suddenly collapses. This process can be responsible for the periodic sand bursts observed in experiments. The model was applied to 2D plain strain geometry for a long wellbore or perforation using a finite difference software. The 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 flow direction. The resistance forces are the forces that retain the elements in place and are generated by vertical and tangential stresses and the friction coefficient. The advantage of using this method is that no calibration parameter for sanding (such as sand production coefficient) is necessary. The 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 conflicting assumptions. Material mass balance equations are based on rigid porous media while equilibrium equations are based on deformable porous media. Therefore 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 fluid flow. Therefore, 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.


#### 2.1.3. Phases Involved in Continuum Approaches

The models can be categorized into two groups based on the phases involved. In the first group, mass balance equation is solved for only fluid and solid phases while the second group recognizes fluidized solid as a phase and solves for solid concentration in the fluid. Fluidized solids are the particles in suspension that move with the fluid. 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 fluid phases) and couple the equations with an erosion model. This is mainly because good agreement between the model and field observations were obtained when combined with a suitable sanding criterion.Multiphase fluid flow 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 inflow 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 flow (oil and water) and 3D geometry. In this work, the effect of water on rock strength reduction is reflected in material properties such as cohesion. The results show that water contact has a significant impact on the sand rate.Gas flow also may hasten the instability process in sand production. When gas comes out the oil phase due to pressure drop and flows 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 finite element method. They assumed that eroded mass in erosion law is proportional to the total fluid flux, which is referred to the oil and gas flux. The effect of multiphase fluid flow was also considered by Wang and Xue [32] for oil-water phases and later by Wang et al. [35] for oil-water-gas phases.The 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. The first one is to remove such elements from the model assuming cavities and wormholes grow as a result of sanding. This seems to be a suitable approach in stronger rocks in which stable cavities can form. The 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 infill material, or cohesionless sand particles. In this approach, sandstone properties are altered to those of degraded cohesionless sand. The 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. There 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. This is mainly attributed to the sand deposition and plugging of the pore space, which is not the case for less permeable sand [26].The 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). These numbers are the lowest values reported for loose sand [70].


#### 2.1.5. Model Design

Openhole completion is often 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,σh and σH, 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 finally higher sanding rate. One solution is using 3D simulation but it is computationally demanding because very fine 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 filled with infill 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 influences the pressure gradient around the perforations and also impacts the mechanical response. These models are also unable to capture the perforation direction effect, which can play a significant role in perforation stability and failure as demonstrated by Papanastasiou and Zervos [71]. They 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. These models require calibration against field and/or physical model testing before application to real-world sanding problems.


#### 2.1.6. 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].The current models are unable of predicting the sand bridging and fines retention in the rock. Sand particles can aggregate at the perforation cavity and form a stable entity calledsand bridge and act as a filter which may reduce further sand production as long as the flow rate remains constant. The 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. The difficult part about modeling sand deposition is the calibration of the critical porosity or the critical sand concentration after 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. The observation in injection wells is often described as the cases with high sand production within a short duration. Few papers [72–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 defined 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 flow easily like a liquid. As no investigation has been performed on liquefaction in sand production, it is difficult to confirm its role in massive sand production.


### 2.2. 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. Therefore, discontinuum approach is promising to simulate phenomena such as detachment of individual particles from the rock matrix.Cundall [83] first introduced the Discrete Element Method (DEM). The 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. The contact forces and displacements of a stressed assembly of particles are found by tracing the movements of the individual particles [84].Table3 summarizes some of the discontinuum-based sanding models. At first, O’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. They used particle bonding scheme to mimic the cement and cohesion due to capillary forces. The bond also incorporated spring stiffness and a nominal tensile breaking strength assuming a bond dimension proportional to the size of the particles it connects. They incorporated the fluid flow calculations by combining the continuity equation and Darcy’s law using the finite element method. Darcy flow is formulated with a measure for effective permeability in the solid medium based on the porosity and average diameter of the solid particles. Their 2D model provides an understanding of the fundamental physics involved in sand production and the relative importance of various rock and fluid properties.Table 3
Summary of the numerical works on sand production in the literature (Discontinuum approach).
Model
Geometry and software
Particle shape
Failure criterion of the particle
Fluid flow analysis
O’Connor et al. [77]
2D
Irregular particles
Tensile failure
2D, FE, Darcy’s law
Jensen and Preece [78]
2D, MIMES software
n-sided polygon
Tensile failure
2D, FE, Darcy’s law
Li et al. [79]
2D, PFC2D software
Circular
(i) Tensile failure(ii) Shear failure(iii) Compressive failure
2D, Explicit FD, Darcy’s law
Li and Holt [80]
2D, PFC2D software
Circular
(i) Tensile failure(ii) Shear failure
2D, fluid flow networks, Darcy’s law
Cheung [81]
3D, PFC3D software
Spherical
(i) Tensile failure(ii) Shear failure
1D, Navier-Stokes equations
Zhou et al. [82]
3D, PFC3D software
Spherical
(i) Tensile failure(ii) Shear failure
3D, Navier-Stokes equations (CFD software)Jensen and Preece [78] explored the coupling of 2D DEM and finite element implementation of the 2D continuity equation for Darcy flow to assess the sanding potential. The basic particle shape used by the model was an n-sided polygon and only the tensile mode for bond failure was considered. They 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 fluid flow to study sanding. PFC2D simulates an assembly of circular disks with the bonds inserted between them. The disks are assumed to be rigid but they can overlap. The 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. The 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. They used a failure criterion composed of tensile failure, shear failure, and compressive failure. A simple approach was used to calculate and couple the fluid flow. They found three typical failure patterns in the simulations similar to those observed in laboratory experiments. The 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 fluid flow in 2D DEM codes by introducing fluid flow networks and simulating flow along the flow paths connecting the voids which are referred to as pipes. The fluid velocities that flow through the pipes and the pore pressures were computed based on Darcy’s theory. The forces arising from the pore fluid were then calculated and applied to the particles in the DEM model. Li and Holt [80] implemented this type of fluid-solid coupling system in the PFC2D codes. While the geometrical limitations of using a 2D model to investigate breakout geometry are obvious, the fluid flow simulation through flow networks is computationally expensive.The entire coupling techniques described above used Darcy’s law to calculate fluid flux or pressure. Darcy’s law has been derived from the Navier-Stokes equations via homogenization that is only valid for slow and viscous flow. The conditions may not be met around the wellbore where breakout forms [81]. Chan and Tipthavonnukul [86] proposed a method to couple continuum and discontinuum flow to simulate granular particles movement in a flowing fluid. The fluid flow is modeled using the 2D Navier-Stokes equations solved by a finite volume method. The coupling is achieved by detecting the presence of the solid in the flow domain and altering the flow resistance accordingly. The 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. The 2D DEM models overestimate the effect of fluid flow on the integrity of the assembly of particles as the resistance to particle dislodgment due to contact forces and bonds normal to the fluid flow is neglected in 2D models. In addition, 2D models cannot represent three-dimensional pore flow networks.Cheung [81] used a coupled fluid-solid 3D DEM model for a perforation test simulation to study the sand production problem. They used the Navier-stokes equations assuming radial flow. Although this scheme is computationally inexpensive, it has two major problems. First, by assuming radial flow, it is not suitable for investigating the impact of the flow at the tip of the perforation where the fluid flow is in all directions. Second, the fluid flow scheme does not consider the presence of the cement between particles. The magnitude of the radial pore fluid velocity in each fluid cell is calculated, considering only the presence of the particles. The presence of cement can highly affect the rock conductivity and therefore the fluid velocity.Later, Zhou et al. [82] employed DEM with computational fluid dynamics (CFD) and showed that the main features of sand erosion can be captured by the CFD-DEM approach.The 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. There may be several variations in the parameters and it is difficult to conclude which set of parameters is most appropriate for the material.


### 2.3. 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. Continuum-based approaches can be used in far field 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 field-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 fluid 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.


## 2.1. 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.Table1 summarizes the majority of continuum-based sanding models. Initially many researchers first 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 filtration theory without solving the equilibrium equation. Later Papamichos and Stavropoulou [8] combined the evolution of localized deformation with hydrodynamic erosion. This was the beginning for many researches that adopt full strength hardening/softening behavior of sandstone in their models [9–17]. The results are highly mesh-dependent for strain softening 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 micro-structure 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.Table 1
Summary of the numerical works on sand production in the literature (continuum approach).
Model
Geometry and solution method
Yield
Hardening/softening
Coupling
Phases
Sanding criteria
Permeability alteration in the sanded zone
Other features
Morita et al. [25], Morita et al. [26], Burton et al. [27]
3D; finite element (FE)(SAND3D software)
Kinematic model with a cap
yes (flow friction)
Iteratively  coupled
(1) Fluid(2) Solid
Maximum plastic strain limit
No change
(1) Only the onset of sanding(2) Burton applied it for gas reservoirs
(Vardoulakis et al. [7])
1D; finite difference (FD)
N/A
N/A
Fluid flow and erosion are coupled
(1) Fluid(2) Fluidized    solid
Erosion
k
/
k
0
=
φ
3
/
(
1
-
φ
2
)  (Carman-Kozeny)
(1) Only hydromechanical effects; equilibrium eqn. is not solved(2) Sand deposition is neglected in modeling
Skjaerstein et al. [28]
1D; FD
N/A
N/A
Fluid flow and erosion are coupled.
(1) Fluid(2) Fluidized solid
Erosion
1
/
k
=
B
φ
+
A  (found experimentally)
Forcheimer’s law was used instead of Darcy’s law to account for turbulence
Papamichos and Malmanger  [29]
2D Axial symmetry and 3D; FE;  Newoton-Paphson (NR) iterations
MC
Yes
Fully coupled
(1) Fluid(2) Solid(3) Fluidized solid
Erosion
k
/
k
0
=
φ
3
/
ρ
-
(
1
-
φ
2
)  ρ-=ρf+c(ρs-ρf)
Tension cut-off: function of both plastic strain and porosity; by the factor(1-φ)/(1-φi)
Papamichos et al. [9]
2D Axial symmetry;  FEM; NR iterations
MC
Yes
Fully coupled
(1) fluid(2) solid
Erosion
Carman-Kozeny
Tension cut-off and Young mod changed by the factor(1-φ)/(1-φi)
Yi [30]
2D axial symmetry;  FD
MC
No
(1) Fluid(2) Solid(3) Fluidized sand
Erosion
k
/
k
0
=
γ
(
φ
3
(
1
-
φ
0
2
)
/
φ
0
3
(
1
-
φ
2
)
)
×
(
1
/
(
1
+
β
m
d
/
ρ
s
)
)
Sand deposition in porous media is considered
Vaziri et al. [31]
2D axial symmetry; FE; fully-implicit
Modified MC with tensile failure
Yes
Fully coupled
(1) Fluid(2) Solid
Tensile failure
k
/
k
0
=
exp
[
-
2.88
×
10
-
3
(
σ
i
′
-
σ
0
′
)
]
Zero stiffness, compressibility and highk for the liquefied tensile zone
Wang and Xue [32]
FE; Crank-Nicholson for time integration
MC
No
Fully coupled
(1) Oil(2) Water(3) Solid(4) Fluidized sand
Erosion
Kozeny-Poiseuille law and Carman-Kozeny
Chin and Ramos [33]
2D and 3D; FE; explicit; NR iteration
Drucker-Prager
No
Coupled
(1) Fluid(2) Solid
Shear dilation
Power law with porosity (exponent = 5.6)
Porosity is changed as a function of plastic volumetric strain
Nouri et al. [34]
2D plane strain;  FD
Bilinear MC
Yes
Fully coupled
(1) Fluid(2) Solid
Tensile failure or shear-failed element falls in tension
0
Capillary is considered as a residual cohesion
Wang et al. [35]
2D; FE
Drucker-Prager
No
Fully coupled
(1) Solid(2) Fluidized solid  (3) Oil(4) Water(5) Gas
k
=
k
0
exp
[
A
(
(
φ
-
φ
0
)
/
(
φ
max
-
φ
0
)
)
]
Cohesion and friction drop linearly with porosity
Servant et al. [36]
2D; FE
MC
No
Iteratively  coupled
(1) Fluid(2) Solid(3) Fluidized solid
Yielding
Not mentioned
Failed material is treated as a Poiseuille fluid.Constant viscosity for the slurry
Detournay et al. [11]
2D; FD
MC
Yes
Iteratively  coupled
(1) Fluid(2) Solid
Erosion
k
/
k
0
=
φ
3
(
1
-
φ
0
2
)
/
φ
0
3
(
1
-
φ
2
)
Bulk mod. Change byφ0/φ
Nouri et al. [12]
2D axial symmetry; FD
Bilinear MC
Yes
Iteratively  coupled
(1) Fluid(2) Solid
Tensile failure or a shear-failed element falls in tension
k
=
C
φ
3
/
γ
f
(
1
-
φ
2
)
Nouri et al. [37]
FE; NR iterations
MC
Yes
Fully coupled
(1) Fluid(2) Solid
Tension
0
Adaptive mesh is used
Nouri et al. [38], Vaziri et al. [14]
2D plane strain and axial symmetry; FD
Bilinear MC
Yes
Iteratively  coupled
(1) Fluid(2) Solid
Complete degradation and tensile mean effective stress
High permeability is assigned to infill materials (elements that satisfy sanding criteria)
WH pulsing is included in the model, and stiffness changes with sanding (Vaziri et al. [14])
Detournay [15]
2D plane strain;  FD
Double-yield
Yes
Iteratively  coupled
(1) Fluid(2) Solid
Same as (Detournay et al. [11])
Not mentioned
Only cap yielding is considered
Nouri et al. [16]
All the features are the same as (Nouri et al. [38]) but fracture energy regularization is applied to remove mesh dependency and the calibration of hardening/softening behavior is discussed thoroughly
Kim et al. [17]
3D;  FD
MC
Yes
Not mentioned
(1) Fluid(2) Solid
Calculated from force balance on the element
N/A (element removal)
No calibration parameter to match the experiment is used and yet a good match is observed
Wang et al. [23]
2D;  FE(ELFEN software)
Soft rock model
Yes
Fully coupled
(1) Fluid(2) Solid
Displacement  (the materials passing the wellface are removed)
K changes as a function of dilation or volumetric strain: k=exp(Aεv)
(1) Adaptive mesh(2) Arbitrary Largrangian Eulerian formulation(3) Fracture Energy Regularization
Azadbakht et al.[39]
2D, FE for fluid and FD for solid
Bilinear MC
Yes
Sequentially coupled
(1) Fluid(2) Solid
Erosion
K changes linearly with porosity
Water effect is considered by reducing the cohesionBulk moduli change linearly with porosityAfter critical porosity, residual properties are assigned to the element as infill material


### 2.1.1. 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. This 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. This model may be used as a quick estimate of sanding onset in relation to production parameters and in situ stresses. The 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]; Peden et al. [47]; Papamichos and Vardoulakis [9]; Detournay et al. [11]; Wan and Wang [48]; Servant et al. [36]; Wang et al. [35]; Wan et al. [49]; Detournay [15]; Wan and Wang [50]; Vaziri et al. [14]; Nouri et al. [13]; Nouri et al. [51]; Nouri et al. [34]; Rahmati et al. [24]; Azadbakht et al. [39]). These 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. Therefore, employment of such models in simulation of material degradation leads to inability in recovering size and dependency of the results on numerical mesh design. This 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. These models also proved efficient in capturing the size effect often observed in laboratory tests performed on thick-walled cylinder samples [53].Vardoulakis et al. [54] used bifurcation theory to predict the failure mode and the critical value of the stress at infinity. They showed that failure depends on stress path and boundary conditions. Tronvoll et al. [55] performed finite 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.The 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 field problems.As shown in Table1, 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] modified the MC model using a bilinear yield function to differentiate sand behavior under low and high confining stresses. The model was later used by Nouri et al. [12], Nouri et al. [38], and Vaziri et al. [14]. The 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. Therefore, Detournay [15] extended Detournay  et al. [11] work by accounting for the compaction mode of failure. The model used Double Yield cap constitutive law to capture compaction bands. The 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.


### 2.1.2. Sanding Criteria Used in Continuum Approaches

Several mechanisms are recognized as responsible for sand production. They are mainly based on shear and tensile failure, critical pressure gradient, critical drawdown pressure, critical plastic strain, and erosion criteria. Table2 summarizes the main sanding criteria used in sand models.Table 2
Sand criteria used in continuum-based sanding models.
Sanding criterion
Formulation
Further developments
Assumptions
Vardoulakis et al. [7]
d
φ
/
d
t
=
m
˙
/
ρ
s
=
λ
(
1
-
φ
)
c
∥
q
i
∥  And with sand deposition:= λ(1-φ)(c-c2/ccr)∥qi/q0∥n
Sand deposition is also considered in some models [28]
It is assumed that the sand in place is fully degraded from the beginning and the production is only due to the hydrodynamic forces. Equilibrium eqn. for solid phase is often ignoredTo initiate the process a very small solid concentration is given as a boundary condition. The results are insensitive to this value as long as it is small
Skjaerstein et al. [28]
d
φ
/
d
t
=
m
˙
/
ρ
s
=
λ
c
(
φ
-
φ
2
/
φ
cr
)
q
x
Papamichos and Malmanger [29]
m
˙
/
ρ
s
=
λ
(
1
-
φ
)
c
q
i
q
i  λ=λ(εp) if εp≥εppeak
The critical flow rate is ignoredSoftening is not simulated
Yi [30]
∂
φ
/
∂
t
=
η
λ
(
1
-
φ
-
(
m
d
/
ρ
s
)
)
v
2
-
λ
d
v
c
s  md˙/ρs=λdvcs+(1-η)(m˙/ρs)
Papamichos et al. [9]
ε
p
>
ε
p
cr  m˙/ρs=λ(1-φ)∥qi∥if φ<φcr, else pore collapse or sudden element removal.λ=λ(εp)
Later models definedεpcr as a function of stresses (Wang and Xue [32]). Critical plastic strain εpcr  was later linked to the state of zero cohesion (Azadbakht et al. [39])
Porosity is assumed to increase until it reaches unity. Usually, it is assumed that eroded elements collapse once their porosities reach a critical level less than unity (say 0.5)
(1) Erosioncriterion
Chin and Ramos [33]
d
φ
/
d
t
=
(
1
-
φ
)
∇
·
v
s
Wang [32, 35]
d
φ
/
d
t
=
m
˙
/
ρ
s
=
λ
(
1
-
φ
)
c
q
i
q
i
Wang et al. [35] assumed λ=λ(εp)
The resisting effect of cohesion was ignored resulting in conservative sanding assessmentSoftening is not considered
Servant et al. [36]
After failure, new material properties are assigned to the failed material (such as zero cohesion). The infill is assumed to behave as a Poiseuille fluid with a viscosity depending on the sand concentration
Detournay et al. [11], Detournay [15]
q
>
q
cr  qcr=qcr(k,φ,Rp,∅,C)  if φ<φcr,∂φ/∂t=λ(1-φ)(q-qcr)(SL/VL)  else element removal
Azadbakht et al. [39]
ε
p
≥
ε
p
cr
and
q
>
q
cr  if φ<φcr  m˙/ρs=λ(1-φ)(qi-qcr)λ=λ(εp)qcr=qcr(k,φ,Rp,∅)else pore collapse or element removal
Some models consider sand production from elements whose cohesion is already reduced but not equal to zero yet (Detournay et al. [11]): qcr=qcr(k,φ,Rp,C,∅)
(2) Tensilecriterion
Vaziri et al. [10], Nouri et al. [12],  Nouri et al. [13], Vaziri et al. [14]
Complete degradation (c=0) and tensile failure
Tensile failure has been replaced by either tensile effective stress or tensile mean effective stress, assuming cohesionless sand has no tensile strengthA more comprehensive one is: shear-failed element with zero cohesion falls in tension or sand fails in tension (Nouri et al. [12])
The element that satisfies sanding criteria is removed from the mesh (assuming cavity can grow adjacent to the wellbore)
(3) Force balanceon elements
Kim et al. [17]
A
r
p
∣
r
=
L
-
A
r
∣
r
=
0
p
w
>
A
r
σ
r
∣
r
=
L
r
=
0
-
A
θ
σ
θ
r
∣
θ
=
0
θ
=
L
-
A
z
σ
z
r
∣
z
=
0
z
=
L
+
μ
f
(
A
θ
σ
θ
∣
θ
=
0
+
A
θ
σ
θ
∣
θ
=
L
+
A
z
σ
z
∣
z
=
0
+
A
z
σ
z
∣
z
=
L
)
For illustration purposes the boundary is shown from 0 to L in all directions.This can be considered as an elaboration of tensile criterion considering the effect of friction coefficient
The element is removed suddenly without allowing the porosity to grow graduallyWhen the effective minimum principal stress is equal to the tensile strength of the formation rock, tensile failure may occur. This 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 after 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. This 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. This 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 flow scouring the surface and normal flow 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 flowing flow 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 flow, axial flow parallel to the perforation 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 flow conditions. In the governing equations, they included Brinkman’s extension of Darcy’s law, which accounts for a smooth transition between channel flow and Darcian flow. The 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 fluid system with the erosion behavior of the solids due to fluid flow. 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. The model is based on nonlinear elastoplasticity, nonlinear stress dependent elasticity, friction hardening and cohesion softening, and single-phase flow fully coupled with geomechanics.Wang et al. [35] also performed a fully coupled single-phase flow analysis based on erosion mechanics using the FEM. They 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 flow rate. They modified the erosion law used by Vardoulakis et al. [7] with the addition of a critical flow rate which is a function of the grain size. They also assumed that sand continuously is produced until a critical porosity is reached beyond which the material suddenly collapses. This process can be responsible for the periodic sand bursts observed in experiments. The model was applied to 2D plain strain geometry for a long wellbore or perforation using a finite difference software. The 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 flow direction. The resistance forces are the forces that retain the elements in place and are generated by vertical and tangential stresses and the friction coefficient. The advantage of using this method is that no calibration parameter for sanding (such as sand production coefficient) is necessary. The 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 conflicting assumptions. Material mass balance equations are based on rigid porous media while equilibrium equations are based on deformable porous media. Therefore 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 fluid flow. Therefore, 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.


### 2.1.3. Phases Involved in Continuum Approaches

The models can be categorized into two groups based on the phases involved. In the first group, mass balance equation is solved for only fluid and solid phases while the second group recognizes fluidized solid as a phase and solves for solid concentration in the fluid. Fluidized solids are the particles in suspension that move with the fluid. 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 fluid phases) and couple the equations with an erosion model. This is mainly because good agreement between the model and field observations were obtained when combined with a suitable sanding criterion.Multiphase fluid flow 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 inflow 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 flow (oil and water) and 3D geometry. In this work, the effect of water on rock strength reduction is reflected in material properties such as cohesion. The results show that water contact has a significant impact on the sand rate.Gas flow also may hasten the instability process in sand production. When gas comes out the oil phase due to pressure drop and flows 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 finite element method. They assumed that eroded mass in erosion law is proportional to the total fluid flux, which is referred to the oil and gas flux. The effect of multiphase fluid flow was also considered by Wang and Xue [32] for oil-water phases and later by Wang et al. [35] for oil-water-gas phases.The 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. The first one is to remove such elements from the model assuming cavities and wormholes grow as a result of sanding. This seems to be a suitable approach in stronger rocks in which stable cavities can form. The 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 infill material, or cohesionless sand particles. In this approach, sandstone properties are altered to those of degraded cohesionless sand. The 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. There 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. This is mainly attributed to the sand deposition and plugging of the pore space, which is not the case for less permeable sand [26].The 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). These numbers are the lowest values reported for loose sand [70].


### 2.1.5. Model Design

Openhole completion is often 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,σh and σH, 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 finally higher sanding rate. One solution is using 3D simulation but it is computationally demanding because very fine 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 filled with infill 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 influences the pressure gradient around the perforations and also impacts the mechanical response. These models are also unable to capture the perforation direction effect, which can play a significant role in perforation stability and failure as demonstrated by Papanastasiou and Zervos [71]. They 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. These models require calibration against field and/or physical model testing before application to real-world sanding problems.


### 2.1.6. 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].The current models are unable of predicting the sand bridging and fines retention in the rock. Sand particles can aggregate at the perforation cavity and form a stable entity calledsand bridge and act as a filter which may reduce further sand production as long as the flow rate remains constant. The 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. The difficult part about modeling sand deposition is the calibration of the critical porosity or the critical sand concentration after 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. The observation in injection wells is often described as the cases with high sand production within a short duration. Few papers [72–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 defined 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 flow easily like a liquid. As no investigation has been performed on liquefaction in sand production, it is difficult to confirm its role in massive sand production.


## 2.1.1. 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. This 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. This model may be used as a quick estimate of sanding onset in relation to production parameters and in situ stresses. The 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]; Peden et al. [47]; Papamichos and Vardoulakis [9]; Detournay et al. [11]; Wan and Wang [48]; Servant et al. [36]; Wang et al. [35]; Wan et al. [49]; Detournay [15]; Wan and Wang [50]; Vaziri et al. [14]; Nouri et al. [13]; Nouri et al. [51]; Nouri et al. [34]; Rahmati et al. [24]; Azadbakht et al. [39]). These 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. Therefore, employment of such models in simulation of material degradation leads to inability in recovering size and dependency of the results on numerical mesh design. This 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. These models also proved efficient in capturing the size effect often observed in laboratory tests performed on thick-walled cylinder samples [53].Vardoulakis et al. [54] used bifurcation theory to predict the failure mode and the critical value of the stress at infinity. They showed that failure depends on stress path and boundary conditions. Tronvoll et al. [55] performed finite 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.The 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 field problems.As shown in Table1, 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] modified the MC model using a bilinear yield function to differentiate sand behavior under low and high confining stresses. The model was later used by Nouri et al. [12], Nouri et al. [38], and Vaziri et al. [14]. The 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. Therefore, Detournay [15] extended Detournay  et al. [11] work by accounting for the compaction mode of failure. The model used Double Yield cap constitutive law to capture compaction bands. The 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.


## 2.1.2. Sanding Criteria Used in Continuum Approaches

Several mechanisms are recognized as responsible for sand production. They are mainly based on shear and tensile failure, critical pressure gradient, critical drawdown pressure, critical plastic strain, and erosion criteria. Table2 summarizes the main sanding criteria used in sand models.Table 2
Sand criteria used in continuum-based sanding models.
Sanding criterion
Formulation
Further developments
Assumptions
Vardoulakis et al. [7]
d
φ
/
d
t
=
m
˙
/
ρ
s
=
λ
(
1
-
φ
)
c
∥
q
i
∥  And with sand deposition:= λ(1-φ)(c-c2/ccr)∥qi/q0∥n
Sand deposition is also considered in some models [28]
It is assumed that the sand in place is fully degraded from the beginning and the production is only due to the hydrodynamic forces. Equilibrium eqn. for solid phase is often ignoredTo initiate the process a very small solid concentration is given as a boundary condition. The results are insensitive to this value as long as it is small
Skjaerstein et al. [28]
d
φ
/
d
t
=
m
˙
/
ρ
s
=
λ
c
(
φ
-
φ
2
/
φ
cr
)
q
x
Papamichos and Malmanger [29]
m
˙
/
ρ
s
=
λ
(
1
-
φ
)
c
q
i
q
i  λ=λ(εp) if εp≥εppeak
The critical flow rate is ignoredSoftening is not simulated
Yi [30]
∂
φ
/
∂
t
=
η
λ
(
1
-
φ
-
(
m
d
/
ρ
s
)
)
v
2
-
λ
d
v
c
s  md˙/ρs=λdvcs+(1-η)(m˙/ρs)
Papamichos et al. [9]
ε
p
>
ε
p
cr  m˙/ρs=λ(1-φ)∥qi∥if φ<φcr, else pore collapse or sudden element removal.λ=λ(εp)
Later models definedεpcr as a function of stresses (Wang and Xue [32]). Critical plastic strain εpcr  was later linked to the state of zero cohesion (Azadbakht et al. [39])
Porosity is assumed to increase until it reaches unity. Usually, it is assumed that eroded elements collapse once their porosities reach a critical level less than unity (say 0.5)
(1) Erosioncriterion
Chin and Ramos [33]
d
φ
/
d
t
=
(
1
-
φ
)
∇
·
v
s
Wang [32, 35]
d
φ
/
d
t
=
m
˙
/
ρ
s
=
λ
(
1
-
φ
)
c
q
i
q
i
Wang et al. [35] assumed λ=λ(εp)
The resisting effect of cohesion was ignored resulting in conservative sanding assessmentSoftening is not considered
Servant et al. [36]
After failure, new material properties are assigned to the failed material (such as zero cohesion). The infill is assumed to behave as a Poiseuille fluid with a viscosity depending on the sand concentration
Detournay et al. [11], Detournay [15]
q
>
q
cr  qcr=qcr(k,φ,Rp,∅,C)  if φ<φcr,∂φ/∂t=λ(1-φ)(q-qcr)(SL/VL)  else element removal
Azadbakht et al. [39]
ε
p
≥
ε
p
cr
and
q
>
q
cr  if φ<φcr  m˙/ρs=λ(1-φ)(qi-qcr)λ=λ(εp)qcr=qcr(k,φ,Rp,∅)else pore collapse or element removal
Some models consider sand production from elements whose cohesion is already reduced but not equal to zero yet (Detournay et al. [11]): qcr=qcr(k,φ,Rp,C,∅)
(2) Tensilecriterion
Vaziri et al. [10], Nouri et al. [12],  Nouri et al. [13], Vaziri et al. [14]
Complete degradation (c=0) and tensile failure
Tensile failure has been replaced by either tensile effective stress or tensile mean effective stress, assuming cohesionless sand has no tensile strengthA more comprehensive one is: shear-failed element with zero cohesion falls in tension or sand fails in tension (Nouri et al. [12])
The element that satisfies sanding criteria is removed from the mesh (assuming cavity can grow adjacent to the wellbore)
(3) Force balanceon elements
Kim et al. [17]
A
r
p
∣
r
=
L
-
A
r
∣
r
=
0
p
w
>
A
r
σ
r
∣
r
=
L
r
=
0
-
A
θ
σ
θ
r
∣
θ
=
0
θ
=
L
-
A
z
σ
z
r
∣
z
=
0
z
=
L
+
μ
f
(
A
θ
σ
θ
∣
θ
=
0
+
A
θ
σ
θ
∣
θ
=
L
+
A
z
σ
z
∣
z
=
0
+
A
z
σ
z
∣
z
=
L
)
For illustration purposes the boundary is shown from 0 to L in all directions.This can be considered as an elaboration of tensile criterion considering the effect of friction coefficient
The element is removed suddenly without allowing the porosity to grow graduallyWhen the effective minimum principal stress is equal to the tensile strength of the formation rock, tensile failure may occur. This 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 after 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. This 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. This 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 flow scouring the surface and normal flow 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 flowing flow 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 flow, axial flow parallel to the perforation 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 flow conditions. In the governing equations, they included Brinkman’s extension of Darcy’s law, which accounts for a smooth transition between channel flow and Darcian flow. The 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 fluid system with the erosion behavior of the solids due to fluid flow. 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. The model is based on nonlinear elastoplasticity, nonlinear stress dependent elasticity, friction hardening and cohesion softening, and single-phase flow fully coupled with geomechanics.Wang et al. [35] also performed a fully coupled single-phase flow analysis based on erosion mechanics using the FEM. They 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 flow rate. They modified the erosion law used by Vardoulakis et al. [7] with the addition of a critical flow rate which is a function of the grain size. They also assumed that sand continuously is produced until a critical porosity is reached beyond which the material suddenly collapses. This process can be responsible for the periodic sand bursts observed in experiments. The model was applied to 2D plain strain geometry for a long wellbore or perforation using a finite difference software. The 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 flow direction. The resistance forces are the forces that retain the elements in place and are generated by vertical and tangential stresses and the friction coefficient. The advantage of using this method is that no calibration parameter for sanding (such as sand production coefficient) is necessary. The 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 conflicting assumptions. Material mass balance equations are based on rigid porous media while equilibrium equations are based on deformable porous media. Therefore 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 fluid flow. Therefore, 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.


## 2.1.3. Phases Involved in Continuum Approaches

The models can be categorized into two groups based on the phases involved. In the first group, mass balance equation is solved for only fluid and solid phases while the second group recognizes fluidized solid as a phase and solves for solid concentration in the fluid. Fluidized solids are the particles in suspension that move with the fluid. 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 fluid phases) and couple the equations with an erosion model. This is mainly because good agreement between the model and field observations were obtained when combined with a suitable sanding criterion.Multiphase fluid flow 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 inflow 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 flow (oil and water) and 3D geometry. In this work, the effect of water on rock strength reduction is reflected in material properties such as cohesion. The results show that water contact has a significant impact on the sand rate.Gas flow also may hasten the instability process in sand production. When gas comes out the oil phase due to pressure drop and flows 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 finite element method. They assumed that eroded mass in erosion law is proportional to the total fluid flux, which is referred to the oil and gas flux. The effect of multiphase fluid flow was also considered by Wang and Xue [32] for oil-water phases and later by Wang et al. [35] for oil-water-gas phases.The 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. The first one is to remove such elements from the model assuming cavities and wormholes grow as a result of sanding. This seems to be a suitable approach in stronger rocks in which stable cavities can form. The 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 infill material, or cohesionless sand particles. In this approach, sandstone properties are altered to those of degraded cohesionless sand. The 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. There 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. This is mainly attributed to the sand deposition and plugging of the pore space, which is not the case for less permeable sand [26].The 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). These numbers are the lowest values reported for loose sand [70].


## 2.1.5. Model Design

Openhole completion is often 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,σh and σH, 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 finally higher sanding rate. One solution is using 3D simulation but it is computationally demanding because very fine 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 filled with infill 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 influences the pressure gradient around the perforations and also impacts the mechanical response. These models are also unable to capture the perforation direction effect, which can play a significant role in perforation stability and failure as demonstrated by Papanastasiou and Zervos [71]. They 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. These models require calibration against field and/or physical model testing before application to real-world sanding problems.


## 2.1.6. 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].The current models are unable of predicting the sand bridging and fines retention in the rock. Sand particles can aggregate at the perforation cavity and form a stable entity calledsand bridge and act as a filter which may reduce further sand production as long as the flow rate remains constant. The 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. The difficult part about modeling sand deposition is the calibration of the critical porosity or the critical sand concentration after 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. The observation in injection wells is often described as the cases with high sand production within a short duration. Few papers [72–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 defined 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 flow easily like a liquid. As no investigation has been performed on liquefaction in sand production, it is difficult to confirm its role in massive sand production.


## 2.2. 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. Therefore, discontinuum approach is promising to simulate phenomena such as detachment of individual particles from the rock matrix.Cundall [83] first introduced the Discrete Element Method (DEM). The 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. The contact forces and displacements of a stressed assembly of particles are found by tracing the movements of the individual particles [84].Table3 summarizes some of the discontinuum-based sanding models. At first, O’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. They used particle bonding scheme to mimic the cement and cohesion due to capillary forces. The bond also incorporated spring stiffness and a nominal tensile breaking strength assuming a bond dimension proportional to the size of the particles it connects. They incorporated the fluid flow calculations by combining the continuity equation and Darcy’s law using the finite element method. Darcy flow is formulated with a measure for effective permeability in the solid medium based on the porosity and average diameter of the solid particles. Their 2D model provides an understanding of the fundamental physics involved in sand production and the relative importance of various rock and fluid properties.Table 3
Summary of the numerical works on sand production in the literature (Discontinuum approach).
Model
Geometry and software
Particle shape
Failure criterion of the particle
Fluid flow analysis
O’Connor et al. [77]
2D
Irregular particles
Tensile failure
2D, FE, Darcy’s law
Jensen and Preece [78]
2D, MIMES software
n-sided polygon
Tensile failure
2D, FE, Darcy’s law
Li et al. [79]
2D, PFC2D software
Circular
(i) Tensile failure(ii) Shear failure(iii) Compressive failure
2D, Explicit FD, Darcy’s law
Li and Holt [80]
2D, PFC2D software
Circular
(i) Tensile failure(ii) Shear failure
2D, fluid flow networks, Darcy’s law
Cheung [81]
3D, PFC3D software
Spherical
(i) Tensile failure(ii) Shear failure
1D, Navier-Stokes equations
Zhou et al. [82]
3D, PFC3D software
Spherical
(i) Tensile failure(ii) Shear failure
3D, Navier-Stokes equations (CFD software)Jensen and Preece [78] explored the coupling of 2D DEM and finite element implementation of the 2D continuity equation for Darcy flow to assess the sanding potential. The basic particle shape used by the model was an n-sided polygon and only the tensile mode for bond failure was considered. They 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 fluid flow to study sanding. PFC2D simulates an assembly of circular disks with the bonds inserted between them. The disks are assumed to be rigid but they can overlap. The 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. The 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. They used a failure criterion composed of tensile failure, shear failure, and compressive failure. A simple approach was used to calculate and couple the fluid flow. They found three typical failure patterns in the simulations similar to those observed in laboratory experiments. The 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 fluid flow in 2D DEM codes by introducing fluid flow networks and simulating flow along the flow paths connecting the voids which are referred to as pipes. The fluid velocities that flow through the pipes and the pore pressures were computed based on Darcy’s theory. The forces arising from the pore fluid were then calculated and applied to the particles in the DEM model. Li and Holt [80] implemented this type of fluid-solid coupling system in the PFC2D codes. While the geometrical limitations of using a 2D model to investigate breakout geometry are obvious, the fluid flow simulation through flow networks is computationally expensive.The entire coupling techniques described above used Darcy’s law to calculate fluid flux or pressure. Darcy’s law has been derived from the Navier-Stokes equations via homogenization that is only valid for slow and viscous flow. The conditions may not be met around the wellbore where breakout forms [81]. Chan and Tipthavonnukul [86] proposed a method to couple continuum and discontinuum flow to simulate granular particles movement in a flowing fluid. The fluid flow is modeled using the 2D Navier-Stokes equations solved by a finite volume method. The coupling is achieved by detecting the presence of the solid in the flow domain and altering the flow resistance accordingly. The 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. The 2D DEM models overestimate the effect of fluid flow on the integrity of the assembly of particles as the resistance to particle dislodgment due to contact forces and bonds normal to the fluid flow is neglected in 2D models. In addition, 2D models cannot represent three-dimensional pore flow networks.Cheung [81] used a coupled fluid-solid 3D DEM model for a perforation test simulation to study the sand production problem. They used the Navier-stokes equations assuming radial flow. Although this scheme is computationally inexpensive, it has two major problems. First, by assuming radial flow, it is not suitable for investigating the impact of the flow at the tip of the perforation where the fluid flow is in all directions. Second, the fluid flow scheme does not consider the presence of the cement between particles. The magnitude of the radial pore fluid velocity in each fluid cell is calculated, considering only the presence of the particles. The presence of cement can highly affect the rock conductivity and therefore the fluid velocity.Later, Zhou et al. [82] employed DEM with computational fluid dynamics (CFD) and showed that the main features of sand erosion can be captured by the CFD-DEM approach.The 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. There may be several variations in the parameters and it is difficult to conclude which set of parameters is most appropriate for the material.


## 2.3. 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. Continuum-based approaches can be used in far field 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 field-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 fluid 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.


## 3. Conclusions

Despite the numerous efforts in sand production and modeling, there are still some fundamental deficiencies which require to be addressed. Considering the works reported in the literature, there is still significant room for improvement in sanding models. Some are listed below.The assessment of the properties of infill materials requires more investigations and experimental data. The choice of these properties plays a significant role in sanding prediction. More accurate correlations for the changes of rock and flow properties with sand production and increase of porosity are essential.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 filter that prevents further sand particles to produce as long as the flow rate remains constant.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 flow, or whether the produced material is sand or shale/clay.The calibration procedure in the DEM model for determining micro material parameters needs further research.In order to have more accurate analysis in the DEM model, the fluid flow scheme needs to be modified. 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.Hybrid model combining DEM for the rocks around the wellbore and continuum approach for far-field rocks is expected to provide a more realistic and yet practical representation of a number of critical factors governing sanding.



---

*Source: 864981-2013-01-09.xml*