# Estimation of Sand Production Rate Using Geomechanical and Hydromechanical Models

Authors: Son Tung Pham
Journal: Advances in Materials Science and Engineering
 (2017)

---

## Abstract

This paper aims to develop a numerical model that can be used in sand control during production phase of an oil and gas well. The model is able to predict not only the onset of sand production using critical bottom hole pressure inferred from geomechanical modelling, but also the mass of sand produced versus time as well as the change of porosity versus space and time using hydromechanical modelling. A detailed workflow of the modelling was presented with each step of calculations. The empirical parameters were calibrated using laboratory data. Then the modelling was applied in a case study of an oilfield in Cuu Long basin. In addition, a sensitivity study of the effect of drawdown pressure was presented in this paper. Moreover, a comparison between results of different hydromechanical models was also addressed. The outcome of this paper demonstrated the possibility of modelling the sand production mass in real cases, opening a new approach in sand control in petroleum industry.


---

## Body

## 1. Introduction

Sand production occurs in many oil fields across the world and it is especially common in the porous sediments. Sand production observed on the surface occurs as a series of three events that happen at downhole area: (1) formation failure, (2) sand erosion due to flow, and (3) sand transport.(1) Formation Failure. In situ stresses and pore pressure act on formation sands and under certain conditions; the criteria for failure are met. The presence of the wellbore and perforations causes a concentration of stresses near these cavities and deformation and failure can occur under certain well-known conditions. This criterion is bottomhole pressure that makes the maximum effective tangential compressive stress equal or higher than the rock strength (failure criteria); significant sanding begins at some point (the onset).(2) Sand Erosion due to Flow. Damaged regions that have failed (meeting the failure criteria) face additional stresses caused by pore pressure gradients. The process of sand erosion is essential for the sand to be removed from the failed region and then to be entrained with the fluid.(3) Sand Transport. Sand erosion detaches sand grains into a perforating cavity or wellbore. Some of the grains are transported to the surface while others settle into the perforation tunnel or into the well hole.During and after sand production, wells can sand-up and that has different effects on the productivity. At first, the productivity seems to increase due to the increase of permeability. However, after a while, the sand produced can obstruct the entrance of hydrocarbon into the wellbore, especially for cased perforated wells using sand screens or gravel pack. Disposal of produced sand is also a significant cost associated with sand production. Finally, sand can be transported to the surface which causes erosion of pipe lines, joints, chokes, and valves. So, if the prediction of sand production is identified, it will help operators to manage the situation properly and prepare suitable treatment methods for the well.This study focuses on the first and second steps of sand production because they have the most impact and they are not well understood. The main reasons why in petroleum industry nowadays we still do not predict the mass of sand production are because of the complexity of numerical models (hence the lack of professional software in this domain) and the unavailable real data of sand production due to the difficulty in collecting this kind of data in the oil field. Therefore, most of the studies predicting the mass of produced sand still stay at the laboratory step. In real life, in petroleum companies’ reports, only Geomechanical models are being used to predict the onset of sand production. This papers aims to bring the application of the Hydromechanical models into a real case in petroleum industry and to combine the use of Geomechanical model, which predicts the onset of sand production, and the Hydromechanical model, which predicts the mass of produced sand.


## 2. Literature Review

Several studies released models predicting the onset of sand production and the amount of sand produced. Parameters affecting sand production have been discussed for decades. However, there is no clear consensus. In this section, a brief review of major conclusions of these past studies is presented.Willson et al. [1] developed a model to predict the sand production rate from the onset of sanding model. The onset of sanding is predicted using a stress–based model of shear failure around a perforation or an open hole. Sand production is assumed to occur once the maximum value of the effective tangential stress around the perforation exceeds the apparent rock strength (rock strength is bf·1.55·TWC strength with bf is boost factor, for cased perforated wells bf = 2; TWC is the Thick Wall Cylinder which is a measure of rock’s strength and is used in sand production study instead of Unconfined Compressive Strength, because TWC reflects more closely the reality of in situ stress sustained by the borehole or the perforation channels than the UCS). No consideration was given to sand transport by drag forces. The model of Willson et al. [1] requires TWC test data for hole size the same as the perforation size. If the perforation is of different size, the application of the model must be careful. Hence the bf boost factor exists to compensate for the difference between the real data and the test data.The criterion for sanding is(1)CBHP≥3σ1-σ2-UCS2-A-poA2-Awhere CBHP is critical bottomhole pressure;σ1 is maximum principal stress; σ2 is minimum principal stress; po is pore pressure. A is poroelastic constant given by A=(1-2v)α/1-v with ν: Poisson ratio and α: Biot’s constant.The model predicts the rate of sand production by utilizing the nondimensionalized concepts of Loading Factor, LF (near-wellbore formation stress normalized by strength), Reynolds number (Re), and water production factor. An empirical relationship between Loading Factor, Reynolds number, and the rate of sand production incorporating the effect of water production was proposed as follows:(2)SPR=fLF,Re,Water-cut;(see [1]).In this formula given by Willson et al. [1], SPR is the sand production rate; water cut is the ratio of water produced compared to the volume of total liquids produced.Although the Willson et al. model [1] takes into account the different phases of the fluid via water cut, the model did not give clear expression of the Sand Production Rate SPR in function of the Loading Factor LF, the Reynolds number, and the water cut.In 1996, Vardoulakis et al. [2] proposed the following sand production model usingmixture theory, assuming that the sand in place is fully degraded from the beginning and the production is due only to the hydrodynamic forces. Equilibrium equation for solid phase is often ignored. The process is initiated with a very small solid concentration given as a boundary condition. The results are insensitive to this value as long as it is small:(3)dϕdt=m˙ρs=λ1-ϕcqfli,where ϕ is the porosity, m˙ is the rate of eroded solid mass per unit volume, ρs is the solid density, λ is the experimentally evaluated sand production coefficient, c is the transport concentration, qfli is the specific discharge in the ith direction, and  is the notation representing the norm of a vector.The Vardoulakis model is difficult to solve because of the complexity of the equations. Moreover, the model does not take into account the different phases of the fluid, so the fluid is considered as single phase, which does not reflect the reality of petroleum fluid. Furthermore, the coefficients were not calibrated due to lack of experimental data.Papamichos et al. [3] developed the following model based on the assumption that failure is due to erosion and porosity increases until it reaches unity. Below is the relation between sand mass and the porosity:(4)m˙ρs=∂ϕ∂t.The dimension of the above equation is thatm˙ is the rate of eroded solid mass per unit volume (g·s−1·m−3), ρs is the solid density (g·m−3), ϕ is the porosity, and t is time in second (s).Variation of sand mass due to erosion is given as(5)m˙ρs=λ1-ϕqfli,where λ is sand production coefficient λ=λ(εp) and  εp is plastic shear strain.The main advantage of Papamichos et al. model [3] is that the authors provided experimental data in order to calibrate the experimental parameters. In this model, not all the plastic deformation areas will produce sand. The sand is produced only when the plastic deformation reaches a limit value (εp>εpeak). However, it is practically difficult to determine this limit value in reality. Besides, the system of differential equations of this model is extremely complex with plenty of empirical parameters.Chin and Ramos [4] developed a sand production model considering that erosionoccurs during sand production as follows, where vs is the solid velocity:(6)dϕdt=1-ϕ∇·vs.The main inconvenient of this model [4] is that it was developed only for weak formation. Besides, the model retains only primary physics of rock failure and coupled rock deformation and fluid flow. The model does not include the effects of well configuration and completion, wellbore storage, erosion, interaction of disaggregated solid and flowing fluid, and solid transport through the porous medium.The analytical model developed by Fjær et al. [5] combined a theoretical model with an empirical relationship as shown below:(7)∂msand∂t=λsandμ1-ϕϕ3qfl-qflcr,where the porosity increases with time:(8)ϕ=ϕ01+4λsandμρsϕ04qfl-qflcrt-t01/4,where qfl is the flow rate, qflcr is the critical flow rate, ϕ0 is the initial porosity, t0 is the initial time, and λsand is proportionality constant; λsand has the dimension of s/m3.Gravanis et al. [6] developed a coupled stress-fluid flow erosion model: (9)m˙ρs=λλpr,t1-ϕβqfl,where λ, with dimension of inverse length, represents the strength of the erosive processes that lead to sand production. β is an exponent parameter of the model. λp(r,t) is an erosion function defined as follows:(10)λpr,t=exp⁡-Λt-aΓ1+1aa×r-rina.Γ ( x ) is the usual Euler Gamma function and Λ is the depth of the plastic region. The profile function λp(r,t) approaches a step function as a increases. The exponent a is fixed at the beginning of the analysis and can be tuned further as part of the calibration procedure. In this work we will fix a=2, which is the same value used in the work of Fjær et al. [5].In 2010, Isehunwa and Olanrewaju [7] proposed a new model for sand production, considering the effects of flow rate, fluid viscosity and density, grain size, and cavity height. Sand is produced by drag and buoyancy forces which predominantly act on the sand particles. The radius of sand production cavity is(11)Ra=qflμ4/9Rs2πHρfg.The volume of sand produced can be expressed as(12)Vsp=πRa2H,where qfl is the fluid flow rate, Rs is the grain radius, H is the cavity height, and g is the gravitational acceleration.Among these models, the ones of Fjær et al. [5] and Gravanis et al. [6] were chosen for this study because of their clarity in the explanation and equations, which are necessary for us to be able to numerically solve the problems. In addition, these models have some advantages in comparison with others. They proposed the basic theory for hydrodynamic erosion of sandstone which is based on filtration theory. They adopted full strength hardening/softening behavior of reservoir stone. They also do take into account the grain size, the gradient elastoplasticity for thick wall cylinders and the cavity failure around boreholes.The models were solved using the workflow developed in Section3 of this paper; then the coding was made using MATLAB software.


## 3. Workflow for Calculation and Calibration

### 3.1. Hydromechanical Erosion Model of Gravanis et al. [6]

According to Gravanis et al. [6], the basic assumptions are as follows:(1)
Fluid flow can be described by Darcy’s law.(2)
We define the mathematical timeT by dT=λ/ΔrAtdt for simple calculations. t is the real time while T is the mathematical time which is introduced to facilitate the resolution of the problem. Δr=rout-rin (Figure 1).
The functionA(t) is related to pressure drop Δp=pout-pin as follows:(13)At=Δpμ∫rinroutdr/kϕr.k is the permeability.(3)
The whole region is divided into a plastic region and an elastic region. In elastic region, we apply Hooke’s law and in plastic region we consider the Mohr-Coulomb failure criterion.(4)
Under the condition that plasticity of the material is damaged and subject to decohesion, it can be eroded under weak hydrodynamic forces. It happens when drawdown pressure (DP) exceeds a critical drawdown pressure (CDP).Figure 1
Schematic of the hollow cylinder [6].


#### 3.1.1. Calculation Workflow

Step 1 (atT=0).
(i)
Calculate functionA(T):(14)AT=Δpμ∫rinroutdr/kϕr, whereA(T) is related to the flow rate qfl(T) by the relation qfl(T)=2πHA(T),  Δp=pout-pin. (ii)
Calculate pressurep(r,T):(15)pr,T=pin+μk0AT∫rinr1-ϕ2ϕ3drr. (iii)
The depth of plastic regionΛ: the elastic and plastic regions are presented in Figure 1 according to Gravanis et al. [6]. In elastic region, combining with the boundary condition σrelastic(rout)=σout, we have(16)σr=σout+C22G1rout2-1r2+α2GλLame+2G1r2∫routrr′pr′dr′σθ=σout+C22G1rout2+1r2+α2GλLame+2Gpr-1r2∫routrr′pr′dr′.r varies from rin to rout. In plastic region, we consider the Mohr-Coulomb failure criterion and combine with the boundary conditionσ(rin)=σin:(17)σr=σinrK-1rinK-1-S0K-11-rK-1rinK-1-αK-1rK-1∫rinrr′-Kpr′dr′σθ=αp1-K+UCS+KσinrK-1rinK-1-S0K-11-rK-1rinK-1-αK-1rK-1∫rinrr′-Kpr′dr′, whereS0 is material cohesion [MPa]. α is Biot’s constant variant from 0 to 1; here we shall set this value to 1 meaning that we will neglect any compressibility effects. K is the principal stress ratio which is equal to σh/σv with σh being the horizontal stress and σv the vertical stress. At this stage, there are two continuity conditions and two unknowns, the integration constantC2 and the location of the plastic zone boundary location R. The conditions are shown below:(18)σrelasticR=σrplasticRσθelasticR=σθplasticR. C 2 is given explicitly in terms of R which is solved numerically for R [6]. Note that (18) are purely algebraic. Thus, the plastic region depth is determined as Λ=R-rin.Step 2 (at timeT=δT).
(i)
The porosity field is calculated from erosion model equation (9):(19)ϕr,T=1-1-ϕ01+β-11-ϕ0β-1Δr/r∫0Tλpr,TdT1/β-1ifβ≠1ϕr,T=1-1-ϕ0exp⁡-Δrr∫0Tλpr,TdTifβ=1. The formula involves the integral∫0Tλpr,TdT which is approximated by λpr,0δT and λp(r,0). Λ0=Λ (T=0) is known from the previous time step where β is an exponent parameter of the model and must be tuned by calibration. (ii)
Calculation of porosity and radius of sand production cavity:(20)ϕ¯T=ΛT-1∫rinroutλpr,Tϕr,Tdr. Radius of sand production cavityra is calculated from wellbore to the distance where porosity of formation is constant and equal to initial porosity ϕ0.(iii)
Calculate functionA(T), pressure p(r,T), and the plastic region depth Λ1 from (14) to (17). (iv)
The erosion strength isλ:(21)λΛ=λ0+λ2-λ0Λ-Λ0Λp-Λ0,Λ≤Λpλ2,Λ≥Λp, whereλ0, λ2, and Λp are model parameters that must be tuned by calibration, λ0 is initial erosion strength [m−1], λ2 is maximum erosion strength [m−1], and Λp is threshold depth of the λ plateau. (v)
Calculate sand production mass rate:(22)M˙=∫m˙dV=ρs·∂ϕ-∂t·π·ra2-rin2·L. Note that this equation must be used in the next time step along withΛp.Step 3 (at the general timeT=iδT).
(i)
The porosity field is calculated from (20) and (21), where the integral ∫0Tλpr,TdT is approximated by λpr,0+⋯+λpr,(i-1)δTδT. This involves the values of the depthΛ at all the previous time steps, which are known. As explained above, the function A, the pressure, and the new value of the plastic region depth Λi=Λ(T=iδΤ) can then be calculated.


### 3.2. Analytical Model of Fjær et al. [5]

The analytical model is based on these assumptions:(1)
The driving mechanism for continuous sand production is erosion from plastified material in the vicinity of the production cavity.(2)
The sand production rate depends on (1) how much the well pressure is reduced below the critical sand production pressure, (2) the fluid flow rate and the fluid viscosity, and (3) the cementation of the rock.(3)
The sand in place is fully degraded from the beginning and the production is due only to the hydrodynamic forces. These forces are proportional to the fluid pressure drop over the volume element, and the pressure drop is proportional to the fluid flow rate as specified by Darcy’s law; moreover the permeability is given by the Kozeny-Carman equation; we have(23)∂msand∂t=λsandμ1-ϕϕ3qfl-qflcr,whereμ is the fluid viscosity. The proportionality constant λsand has the dimension of s/m3(4)
Firstly, it is required that DP > CDP, which expresses the fact that stress induced damage of the rock is a necessary condition for sand production. Secondly, it is required thatqfl>qflcr (critical fluid flow rate), which expresses that the hydrodynamic forces must be strong enough to uphold the erosion process.


#### 3.2.1. Calculation Workflow

Following Fjær et al. [5] we simply assume that the stiffness of the rock remains constant until the porosity has reached a critical value ϕcr. At that point, the entire part of the rock that has been producing sand collapses, and the remaining solid material in that part is produced in one burst. The stresses around the cavity are then redistributed, and the situation is similar to the initial state at t=to, except the cavity has become a little larger.Step 1 (continuous sand production<ϕcr).
(i)
Calculate porosity from (8):(24)ϕ=ϕ01+4λsandμρsϕ04qfl-qflcrt-t01/4. (ii)
Calculate volumetric sand productionVsp: we assume that sand production is related to plastification of the rock, so that only parts of the rock that have suffered some plastic deformations are in a state where sand can be produced.(25)Rcr=RceDP-CDP/UCS≈Rc1+DP-CDPUCSVsp=πRcr2-Rc2L≈2πLRc2DP-CDPUCS, whereRc is radius of the cavity and Rcr is sand production zone (Figure 2). (iii)
Sand production mass:(26)M=ρsϕ-ϕ0VspFigure 2
Illustration of the plastified zone and the sand producing zone around a cavity [8].Step 2 (instantaneous sand productionϕ=ϕcr).
Once the entire part of the rock that has been producing sand collapses, the remaining solid material in that part is produced in one burst.The collapse of the sand producing zone implies that the radius of the cavity increases, fromRc to (1+a1)·Rc, where a1 can be found from (27):(27)a1=RcrRc-1≈DP-CDPUCS(28)Rc⟶1+a1Rc;Rcr⟶1+a1Rcr.(i)
Cumulative sand production: sand mass is given as the initial amount of solid material in the sand producing zone; that is,(29)M=Vspρs1-ϕ0.Now the stresses around the cavity are redistributed, and the situation is the same as it was att=to. Thus, the erosion process starts a new cycle.


### 3.3. Calibration of Empirical Parameters of the Models of Fjær et al. [5] and Gravanis et al. [6]

We use experimental data profile of Papamichos et al. [3] to calibrate empirical parameters (Table 1).Table 1
Experimental data of Papamichos et al. [3].
Variable Value Cylinder internal radius,rin [m] 0.01 Cylinder external radius,rout [m] 0.1 Cylinder height,H [m] 0.2 Internal pore pressure,Pin [Mpa] 0 External pore pressure,Pout [Mpa] 0.15 Inner radial stress,Sin [Mpa] 0 Outer radial stress,Sout [Mpa] 11 Flow rate,qfl [l/min] 0.5 Ratioa1 0.03 Young’s modulus,E [MPa] 6750 Poisson ratio,ν [—] 0.19 Biot’s ratio,α [—] 1 Cohesion,So [MPa] 3.7 Friction angle,Φ  [°] 37.4 Initial porosity,ϕ0 [—] 0.3 Permeability,k [mD] 500 Kozeny-Carman parameter,ko [m2] 8.96E − 12 Solids density,ρsolid [kg/m3] 2640 Dynamic viscosity,μ [MPa·s] 5E09


#### 3.3.1. Calibration of Empirical Parameters of the Gravanis et al. Model

We have empirical parameters:β, λ0, λ2, and Λp and we slightly varied these parameters around the values obtained by trial and error from calibrating with test data of Papamichos et al., 2001 [3]. The calibration results are shown from Figures 3 to 6. For each calibration, the model curve is compared with the experimental curve.Figure 3
Sand mass produced over time while maximum erosion strengthλ2 varies from 0.078 to 0.94 m−1. β=1, λ0=0.42λ2, and Λp=1.1Λ0.From the results presented from Figures3–6, we choose suitable values of empirical parameters which are: the exponent coefficient (β=1), maximum erosion strength (λ2 = 0.088 m−1), initial erosion strength (λ0 = 0.42λ2 m−1) and threshold depth of the λ plateau (Λp=1.1Λo m). The small discrepancy with the experimental data can be neglected.Figure 4
Sand mass produced over time while changing exponent coefficientβ from 0.9 to 1.3. λ2=0.088, λ0=0.42λ2, and Λp=1.1Λ0.Figure 5
Sand mass produced over time while changing ratio of initial to maximum erosion strengthλ0/λ2 from 0.3 to 0.46. λ2=0.088, β=1, and Λp=1.1Λ0.Figure 6
Sand mass produced over time while changing ratio of threshold depth and initial plastic regionΛp/Λ0 from 1.08 to 1.2. β=1, λ2=0.088, and λ0=0.42λ2.


#### 3.3.2. Calibration of Empirical Parameters of Fjær et al. Model

The model of Fjær et al. has three empirical parameters:qflcr, ϕcr, and λsand. Similar to the model of Gravanis et al. [6], we use experimental data of Papamichos et al., 2001 [3] and then run the model and compare with the corresponding experimental curve to calibrate empirical parameters. Results are presented in Figures 7, 8, and 9.Figure 7
Sand mass produced over time while changing critical porosityϕcr from 0.36 to 0.42.Figure 8
Sand mass produced over time while changing critical fluid fluxqflcr from 0.00008 to 0.00014 m/s.Figure 9
Sand mass produced over time coefficient sand productionλsand from 500 to 650 m3/s.Figures7, 8, and 9 indicate that the erosional process gradually increases until critical porosity is reached. Each step is followed by a long period of continuous production at a low rate. When the porosity has reached a critical value ϕcr, at that point, the collapsed zone of the remained solids is produced, and process of sand production is immediate and rapid. Then the erosion process starts on a new cycle.In fact, the result of the model seems very different than the experimental result because Fjaer’s model considers a step of “collapse” (Step2: instantaneous sand production), but overall, the sand production mass over time is matched between the model and the experimental data and that is why we could calibrate the parameters.These results also allow us to choose suitable values of empirical parameters which are critical porosity (ϕcr = 0.4), critical fluid flux (qflcr = 0.001 m/s), and sand production coefficient (λsand = 600 m3/s).


## 3.1. Hydromechanical Erosion Model of Gravanis et al. [6]

According to Gravanis et al. [6], the basic assumptions are as follows:(1)
Fluid flow can be described by Darcy’s law.(2)
We define the mathematical timeT by dT=λ/ΔrAtdt for simple calculations. t is the real time while T is the mathematical time which is introduced to facilitate the resolution of the problem. Δr=rout-rin (Figure 1).
The functionA(t) is related to pressure drop Δp=pout-pin as follows:(13)At=Δpμ∫rinroutdr/kϕr.k is the permeability.(3)
The whole region is divided into a plastic region and an elastic region. In elastic region, we apply Hooke’s law and in plastic region we consider the Mohr-Coulomb failure criterion.(4)
Under the condition that plasticity of the material is damaged and subject to decohesion, it can be eroded under weak hydrodynamic forces. It happens when drawdown pressure (DP) exceeds a critical drawdown pressure (CDP).Figure 1
Schematic of the hollow cylinder [6].


### 3.1.1. Calculation Workflow

Step 1 (atT=0).
(i)
Calculate functionA(T):(14)AT=Δpμ∫rinroutdr/kϕr, whereA(T) is related to the flow rate qfl(T) by the relation qfl(T)=2πHA(T),  Δp=pout-pin. (ii)
Calculate pressurep(r,T):(15)pr,T=pin+μk0AT∫rinr1-ϕ2ϕ3drr. (iii)
The depth of plastic regionΛ: the elastic and plastic regions are presented in Figure 1 according to Gravanis et al. [6]. In elastic region, combining with the boundary condition σrelastic(rout)=σout, we have(16)σr=σout+C22G1rout2-1r2+α2GλLame+2G1r2∫routrr′pr′dr′σθ=σout+C22G1rout2+1r2+α2GλLame+2Gpr-1r2∫routrr′pr′dr′.r varies from rin to rout. In plastic region, we consider the Mohr-Coulomb failure criterion and combine with the boundary conditionσ(rin)=σin:(17)σr=σinrK-1rinK-1-S0K-11-rK-1rinK-1-αK-1rK-1∫rinrr′-Kpr′dr′σθ=αp1-K+UCS+KσinrK-1rinK-1-S0K-11-rK-1rinK-1-αK-1rK-1∫rinrr′-Kpr′dr′, whereS0 is material cohesion [MPa]. α is Biot’s constant variant from 0 to 1; here we shall set this value to 1 meaning that we will neglect any compressibility effects. K is the principal stress ratio which is equal to σh/σv with σh being the horizontal stress and σv the vertical stress. At this stage, there are two continuity conditions and two unknowns, the integration constantC2 and the location of the plastic zone boundary location R. The conditions are shown below:(18)σrelasticR=σrplasticRσθelasticR=σθplasticR. C 2 is given explicitly in terms of R which is solved numerically for R [6]. Note that (18) are purely algebraic. Thus, the plastic region depth is determined as Λ=R-rin.Step 2 (at timeT=δT).
(i)
The porosity field is calculated from erosion model equation (9):(19)ϕr,T=1-1-ϕ01+β-11-ϕ0β-1Δr/r∫0Tλpr,TdT1/β-1ifβ≠1ϕr,T=1-1-ϕ0exp⁡-Δrr∫0Tλpr,TdTifβ=1. The formula involves the integral∫0Tλpr,TdT which is approximated by λpr,0δT and λp(r,0). Λ0=Λ (T=0) is known from the previous time step where β is an exponent parameter of the model and must be tuned by calibration. (ii)
Calculation of porosity and radius of sand production cavity:(20)ϕ¯T=ΛT-1∫rinroutλpr,Tϕr,Tdr. Radius of sand production cavityra is calculated from wellbore to the distance where porosity of formation is constant and equal to initial porosity ϕ0.(iii)
Calculate functionA(T), pressure p(r,T), and the plastic region depth Λ1 from (14) to (17). (iv)
The erosion strength isλ:(21)λΛ=λ0+λ2-λ0Λ-Λ0Λp-Λ0,Λ≤Λpλ2,Λ≥Λp, whereλ0, λ2, and Λp are model parameters that must be tuned by calibration, λ0 is initial erosion strength [m−1], λ2 is maximum erosion strength [m−1], and Λp is threshold depth of the λ plateau. (v)
Calculate sand production mass rate:(22)M˙=∫m˙dV=ρs·∂ϕ-∂t·π·ra2-rin2·L. Note that this equation must be used in the next time step along withΛp.Step 3 (at the general timeT=iδT).
(i)
The porosity field is calculated from (20) and (21), where the integral ∫0Tλpr,TdT is approximated by λpr,0+⋯+λpr,(i-1)δTδT. This involves the values of the depthΛ at all the previous time steps, which are known. As explained above, the function A, the pressure, and the new value of the plastic region depth Λi=Λ(T=iδΤ) can then be calculated.


## 3.1.1. Calculation Workflow

Step 1 (atT=0).
(i)
Calculate functionA(T):(14)AT=Δpμ∫rinroutdr/kϕr, whereA(T) is related to the flow rate qfl(T) by the relation qfl(T)=2πHA(T),  Δp=pout-pin. (ii)
Calculate pressurep(r,T):(15)pr,T=pin+μk0AT∫rinr1-ϕ2ϕ3drr. (iii)
The depth of plastic regionΛ: the elastic and plastic regions are presented in Figure 1 according to Gravanis et al. [6]. In elastic region, combining with the boundary condition σrelastic(rout)=σout, we have(16)σr=σout+C22G1rout2-1r2+α2GλLame+2G1r2∫routrr′pr′dr′σθ=σout+C22G1rout2+1r2+α2GλLame+2Gpr-1r2∫routrr′pr′dr′.r varies from rin to rout. In plastic region, we consider the Mohr-Coulomb failure criterion and combine with the boundary conditionσ(rin)=σin:(17)σr=σinrK-1rinK-1-S0K-11-rK-1rinK-1-αK-1rK-1∫rinrr′-Kpr′dr′σθ=αp1-K+UCS+KσinrK-1rinK-1-S0K-11-rK-1rinK-1-αK-1rK-1∫rinrr′-Kpr′dr′, whereS0 is material cohesion [MPa]. α is Biot’s constant variant from 0 to 1; here we shall set this value to 1 meaning that we will neglect any compressibility effects. K is the principal stress ratio which is equal to σh/σv with σh being the horizontal stress and σv the vertical stress. At this stage, there are two continuity conditions and two unknowns, the integration constantC2 and the location of the plastic zone boundary location R. The conditions are shown below:(18)σrelasticR=σrplasticRσθelasticR=σθplasticR. C 2 is given explicitly in terms of R which is solved numerically for R [6]. Note that (18) are purely algebraic. Thus, the plastic region depth is determined as Λ=R-rin.Step 2 (at timeT=δT).
(i)
The porosity field is calculated from erosion model equation (9):(19)ϕr,T=1-1-ϕ01+β-11-ϕ0β-1Δr/r∫0Tλpr,TdT1/β-1ifβ≠1ϕr,T=1-1-ϕ0exp⁡-Δrr∫0Tλpr,TdTifβ=1. The formula involves the integral∫0Tλpr,TdT which is approximated by λpr,0δT and λp(r,0). Λ0=Λ (T=0) is known from the previous time step where β is an exponent parameter of the model and must be tuned by calibration. (ii)
Calculation of porosity and radius of sand production cavity:(20)ϕ¯T=ΛT-1∫rinroutλpr,Tϕr,Tdr. Radius of sand production cavityra is calculated from wellbore to the distance where porosity of formation is constant and equal to initial porosity ϕ0.(iii)
Calculate functionA(T), pressure p(r,T), and the plastic region depth Λ1 from (14) to (17). (iv)
The erosion strength isλ:(21)λΛ=λ0+λ2-λ0Λ-Λ0Λp-Λ0,Λ≤Λpλ2,Λ≥Λp, whereλ0, λ2, and Λp are model parameters that must be tuned by calibration, λ0 is initial erosion strength [m−1], λ2 is maximum erosion strength [m−1], and Λp is threshold depth of the λ plateau. (v)
Calculate sand production mass rate:(22)M˙=∫m˙dV=ρs·∂ϕ-∂t·π·ra2-rin2·L. Note that this equation must be used in the next time step along withΛp.Step 3 (at the general timeT=iδT).
(i)
The porosity field is calculated from (20) and (21), where the integral ∫0Tλpr,TdT is approximated by λpr,0+⋯+λpr,(i-1)δTδT. This involves the values of the depthΛ at all the previous time steps, which are known. As explained above, the function A, the pressure, and the new value of the plastic region depth Λi=Λ(T=iδΤ) can then be calculated.


## 3.2. Analytical Model of Fjær et al. [5]

The analytical model is based on these assumptions:(1)
The driving mechanism for continuous sand production is erosion from plastified material in the vicinity of the production cavity.(2)
The sand production rate depends on (1) how much the well pressure is reduced below the critical sand production pressure, (2) the fluid flow rate and the fluid viscosity, and (3) the cementation of the rock.(3)
The sand in place is fully degraded from the beginning and the production is due only to the hydrodynamic forces. These forces are proportional to the fluid pressure drop over the volume element, and the pressure drop is proportional to the fluid flow rate as specified by Darcy’s law; moreover the permeability is given by the Kozeny-Carman equation; we have(23)∂msand∂t=λsandμ1-ϕϕ3qfl-qflcr,whereμ is the fluid viscosity. The proportionality constant λsand has the dimension of s/m3(4)
Firstly, it is required that DP > CDP, which expresses the fact that stress induced damage of the rock is a necessary condition for sand production. Secondly, it is required thatqfl>qflcr (critical fluid flow rate), which expresses that the hydrodynamic forces must be strong enough to uphold the erosion process.


### 3.2.1. Calculation Workflow

Following Fjær et al. [5] we simply assume that the stiffness of the rock remains constant until the porosity has reached a critical value ϕcr. At that point, the entire part of the rock that has been producing sand collapses, and the remaining solid material in that part is produced in one burst. The stresses around the cavity are then redistributed, and the situation is similar to the initial state at t=to, except the cavity has become a little larger.Step 1 (continuous sand production<ϕcr).
(i)
Calculate porosity from (8):(24)ϕ=ϕ01+4λsandμρsϕ04qfl-qflcrt-t01/4. (ii)
Calculate volumetric sand productionVsp: we assume that sand production is related to plastification of the rock, so that only parts of the rock that have suffered some plastic deformations are in a state where sand can be produced.(25)Rcr=RceDP-CDP/UCS≈Rc1+DP-CDPUCSVsp=πRcr2-Rc2L≈2πLRc2DP-CDPUCS, whereRc is radius of the cavity and Rcr is sand production zone (Figure 2). (iii)
Sand production mass:(26)M=ρsϕ-ϕ0VspFigure 2
Illustration of the plastified zone and the sand producing zone around a cavity [8].Step 2 (instantaneous sand productionϕ=ϕcr).
Once the entire part of the rock that has been producing sand collapses, the remaining solid material in that part is produced in one burst.The collapse of the sand producing zone implies that the radius of the cavity increases, fromRc to (1+a1)·Rc, where a1 can be found from (27):(27)a1=RcrRc-1≈DP-CDPUCS(28)Rc⟶1+a1Rc;Rcr⟶1+a1Rcr.(i)
Cumulative sand production: sand mass is given as the initial amount of solid material in the sand producing zone; that is,(29)M=Vspρs1-ϕ0.Now the stresses around the cavity are redistributed, and the situation is the same as it was att=to. Thus, the erosion process starts a new cycle.


## 3.2.1. Calculation Workflow

Following Fjær et al. [5] we simply assume that the stiffness of the rock remains constant until the porosity has reached a critical value ϕcr. At that point, the entire part of the rock that has been producing sand collapses, and the remaining solid material in that part is produced in one burst. The stresses around the cavity are then redistributed, and the situation is similar to the initial state at t=to, except the cavity has become a little larger.Step 1 (continuous sand production<ϕcr).
(i)
Calculate porosity from (8):(24)ϕ=ϕ01+4λsandμρsϕ04qfl-qflcrt-t01/4. (ii)
Calculate volumetric sand productionVsp: we assume that sand production is related to plastification of the rock, so that only parts of the rock that have suffered some plastic deformations are in a state where sand can be produced.(25)Rcr=RceDP-CDP/UCS≈Rc1+DP-CDPUCSVsp=πRcr2-Rc2L≈2πLRc2DP-CDPUCS, whereRc is radius of the cavity and Rcr is sand production zone (Figure 2). (iii)
Sand production mass:(26)M=ρsϕ-ϕ0VspFigure 2
Illustration of the plastified zone and the sand producing zone around a cavity [8].Step 2 (instantaneous sand productionϕ=ϕcr).
Once the entire part of the rock that has been producing sand collapses, the remaining solid material in that part is produced in one burst.The collapse of the sand producing zone implies that the radius of the cavity increases, fromRc to (1+a1)·Rc, where a1 can be found from (27):(27)a1=RcrRc-1≈DP-CDPUCS(28)Rc⟶1+a1Rc;Rcr⟶1+a1Rcr.(i)
Cumulative sand production: sand mass is given as the initial amount of solid material in the sand producing zone; that is,(29)M=Vspρs1-ϕ0.Now the stresses around the cavity are redistributed, and the situation is the same as it was att=to. Thus, the erosion process starts a new cycle.


## 3.3. Calibration of Empirical Parameters of the Models of Fjær et al. [5] and Gravanis et al. [6]

We use experimental data profile of Papamichos et al. [3] to calibrate empirical parameters (Table 1).Table 1
Experimental data of Papamichos et al. [3].
Variable Value Cylinder internal radius,rin [m] 0.01 Cylinder external radius,rout [m] 0.1 Cylinder height,H [m] 0.2 Internal pore pressure,Pin [Mpa] 0 External pore pressure,Pout [Mpa] 0.15 Inner radial stress,Sin [Mpa] 0 Outer radial stress,Sout [Mpa] 11 Flow rate,qfl [l/min] 0.5 Ratioa1 0.03 Young’s modulus,E [MPa] 6750 Poisson ratio,ν [—] 0.19 Biot’s ratio,α [—] 1 Cohesion,So [MPa] 3.7 Friction angle,Φ  [°] 37.4 Initial porosity,ϕ0 [—] 0.3 Permeability,k [mD] 500 Kozeny-Carman parameter,ko [m2] 8.96E − 12 Solids density,ρsolid [kg/m3] 2640 Dynamic viscosity,μ [MPa·s] 5E09


### 3.3.1. Calibration of Empirical Parameters of the Gravanis et al. Model

We have empirical parameters:β, λ0, λ2, and Λp and we slightly varied these parameters around the values obtained by trial and error from calibrating with test data of Papamichos et al., 2001 [3]. The calibration results are shown from Figures 3 to 6. For each calibration, the model curve is compared with the experimental curve.Figure 3
Sand mass produced over time while maximum erosion strengthλ2 varies from 0.078 to 0.94 m−1. β=1, λ0=0.42λ2, and Λp=1.1Λ0.From the results presented from Figures3–6, we choose suitable values of empirical parameters which are: the exponent coefficient (β=1), maximum erosion strength (λ2 = 0.088 m−1), initial erosion strength (λ0 = 0.42λ2 m−1) and threshold depth of the λ plateau (Λp=1.1Λo m). The small discrepancy with the experimental data can be neglected.Figure 4
Sand mass produced over time while changing exponent coefficientβ from 0.9 to 1.3. λ2=0.088, λ0=0.42λ2, and Λp=1.1Λ0.Figure 5
Sand mass produced over time while changing ratio of initial to maximum erosion strengthλ0/λ2 from 0.3 to 0.46. λ2=0.088, β=1, and Λp=1.1Λ0.Figure 6
Sand mass produced over time while changing ratio of threshold depth and initial plastic regionΛp/Λ0 from 1.08 to 1.2. β=1, λ2=0.088, and λ0=0.42λ2.


### 3.3.2. Calibration of Empirical Parameters of Fjær et al. Model

The model of Fjær et al. has three empirical parameters:qflcr, ϕcr, and λsand. Similar to the model of Gravanis et al. [6], we use experimental data of Papamichos et al., 2001 [3] and then run the model and compare with the corresponding experimental curve to calibrate empirical parameters. Results are presented in Figures 7, 8, and 9.Figure 7
Sand mass produced over time while changing critical porosityϕcr from 0.36 to 0.42.Figure 8
Sand mass produced over time while changing critical fluid fluxqflcr from 0.00008 to 0.00014 m/s.Figure 9
Sand mass produced over time coefficient sand productionλsand from 500 to 650 m3/s.Figures7, 8, and 9 indicate that the erosional process gradually increases until critical porosity is reached. Each step is followed by a long period of continuous production at a low rate. When the porosity has reached a critical value ϕcr, at that point, the collapsed zone of the remained solids is produced, and process of sand production is immediate and rapid. Then the erosion process starts on a new cycle.In fact, the result of the model seems very different than the experimental result because Fjaer’s model considers a step of “collapse” (Step2: instantaneous sand production), but overall, the sand production mass over time is matched between the model and the experimental data and that is why we could calibrate the parameters.These results also allow us to choose suitable values of empirical parameters which are critical porosity (ϕcr = 0.4), critical fluid flux (qflcr = 0.001 m/s), and sand production coefficient (λsand = 600 m3/s).


## 3.3.1. Calibration of Empirical Parameters of the Gravanis et al. Model

We have empirical parameters:β, λ0, λ2, and Λp and we slightly varied these parameters around the values obtained by trial and error from calibrating with test data of Papamichos et al., 2001 [3]. The calibration results are shown from Figures 3 to 6. For each calibration, the model curve is compared with the experimental curve.Figure 3
Sand mass produced over time while maximum erosion strengthλ2 varies from 0.078 to 0.94 m−1. β=1, λ0=0.42λ2, and Λp=1.1Λ0.From the results presented from Figures3–6, we choose suitable values of empirical parameters which are: the exponent coefficient (β=1), maximum erosion strength (λ2 = 0.088 m−1), initial erosion strength (λ0 = 0.42λ2 m−1) and threshold depth of the λ plateau (Λp=1.1Λo m). The small discrepancy with the experimental data can be neglected.Figure 4
Sand mass produced over time while changing exponent coefficientβ from 0.9 to 1.3. λ2=0.088, λ0=0.42λ2, and Λp=1.1Λ0.Figure 5
Sand mass produced over time while changing ratio of initial to maximum erosion strengthλ0/λ2 from 0.3 to 0.46. λ2=0.088, β=1, and Λp=1.1Λ0.Figure 6
Sand mass produced over time while changing ratio of threshold depth and initial plastic regionΛp/Λ0 from 1.08 to 1.2. β=1, λ2=0.088, and λ0=0.42λ2.


## 3.3.2. Calibration of Empirical Parameters of Fjær et al. Model

The model of Fjær et al. has three empirical parameters:qflcr, ϕcr, and λsand. Similar to the model of Gravanis et al. [6], we use experimental data of Papamichos et al., 2001 [3] and then run the model and compare with the corresponding experimental curve to calibrate empirical parameters. Results are presented in Figures 7, 8, and 9.Figure 7
Sand mass produced over time while changing critical porosityϕcr from 0.36 to 0.42.Figure 8
Sand mass produced over time while changing critical fluid fluxqflcr from 0.00008 to 0.00014 m/s.Figure 9
Sand mass produced over time coefficient sand productionλsand from 500 to 650 m3/s.Figures7, 8, and 9 indicate that the erosional process gradually increases until critical porosity is reached. Each step is followed by a long period of continuous production at a low rate. When the porosity has reached a critical value ϕcr, at that point, the collapsed zone of the remained solids is produced, and process of sand production is immediate and rapid. Then the erosion process starts on a new cycle.In fact, the result of the model seems very different than the experimental result because Fjaer’s model considers a step of “collapse” (Step2: instantaneous sand production), but overall, the sand production mass over time is matched between the model and the experimental data and that is why we could calibrate the parameters.These results also allow us to choose suitable values of empirical parameters which are critical porosity (ϕcr = 0.4), critical fluid flux (qflcr = 0.001 m/s), and sand production coefficient (λsand = 600 m3/s).


## 4. Case Study

An application for the oilfield X in Cuu Long basin using Geomechanical model of Willson et al. [1] to predict the onset of sand production condition is now presented. In addition, it was combined with the Hydromechanical model of Fjær et al. [5] and Gravanis et al. [6] to predict sand production rate; then the results of the models will be compared.After collecting and processing data by analyzing log and other parameters, input data are shown in Table2. The petroleum industry, for historical reasons, has been landed with a mixture of US, British, and SI units, which is often referred to by “oilfield units.” In Vietnam they follow this tradition of using oilfield units, which is a requirement by the government; hence in Table 2 we do not use SI units.Table 2
Data of Well X1.
Input data Value TVD (ft) 9479.5 I (degree) 52 W a z (degree) 244 σ ν (psia) 9300 σ h (psia) 8670 σ H (psia) 9448 ν 0.3 UCS (psia) 2450 α 1 σ h z (degree) 193 ϕ p e r f (degree) 90 TWC (psia) 4936 k (mD) 60 ϕo 0.3 μ (cP) 3 r p e r f = r i n (ft) 0.0328 l p e r f (ft) 1.6405 E (psia) 979004


### 4.1. Calculate Critical Bottomhole Pressure Using Geomechanical Model of Willson et al. [1]

The critical bottomhole pressure CBHP is important information during production phase because it indicates the lowest bottomhole pressure for sand production to not occur. For a specific well in oil field, the only production data that we control is the bottomhole pressure, which is adjusted using surface chokes. The pressure is controlled; hence the flowrate is controlled. For this reason, before predicting the sand production using Hydromechanical models of Gravanis et al. and Fjær et al., we firstly calculate the CBHP using Geomechanical model of Willson et al.The testing showed that a relationship between the effective in situ strength of the formation U and the TWC (Thick wall cylinder) strength relative to a specimen with an OD/ID ratio (Outer/Inner diameters) would be equivalent to(30)UCS=bf∗1.55∗TWC,where bf the boost factor is often taken as bf=2 for cased–hole perforated wells. Equation (1) considers effect of reservoir pressure decline [9] which can be accounted for by updating in situ stresses, while the vertical stresses are usually kept constant:(31)σh=σh+A·ΔPσH=σH+A·ΔP,where A=α1-2ν/1-ν with α is the Biot’s constant and ν is coefficient of Poisson, ΔP=Pc-Pi with Pc is current reservoir pressure, and Pi is initial reservoir pressure.Figure10 presents the bottomhole pressure Pw versus reservoir pressure Pr. In Figure 10 we also present the CBHP when the UCS changes. For example, for UCS = 1950 psia, if reservoir pressure Pr is 6000 psia, the CBHP is about 1500 psia; if the bottomhole pressure Pw is lower than this value 1500 psia, sand production will occur.Figure 10
Prediction of sand production with variation of UCS (Unconfined Compressive Strength).The results showed that if UCS gets smaller, the sand production window is greater. This can be explained by the reason that the formation rocks around the wellbore are weakened when UCS decreases. Thus, risk of sand production increases.We consider a case study with UCS = 2450 psi,Pr = 4000 psi. From Figure 10, we see that if the current reservoir pressure is 4000 psi, then the critical bottomhole pressure is 2020 psia. For sand production to not occur, the bottomhole pressure must be greater than 2020 psi. Otherwise, to increase the flow rate, we must increase DP and if it is greater than CDP, it means that bottomhole pressure must be lower than critical bottomhole pressure (1000 psi); then the sand will be produced. So there is a compromise between increasing flowrate (production requirement) and sand production.


### 4.2. Calculate Sand Production Mass Using the Models of Gravanis et al. and Fjær et al.

We use Hydromechanical Erosion models of Fjær et al. [5] and Gravanis et al. [6] to calculate sand production mass versus time. However, we do not have experimental data of Well X1 to calibrate empirical parameters, so we must use empirical parameters already calibrated in Section 3. In Figure 11 the results of sand mass produced over time calculated by the two models are presented. We also did a sensitivity study of the effect of drawdown pressure. Drawdown pressure, which is the difference between the reservoir pressure and the bottomhole pressure, is the most important data that we control during production phase. We control the drawdown; consequently, we control the production rate. The drawdown is typically controlled by surface chokes. The more the choke is open, the bigger the bottomhole pressure is and, hence, the smaller the drawdown is, and vice versa.Figure 11
Results of two models in the case study of Well X1.Figure11 shows that sand production mass increases when drawdown pressure increases. However, results given by the model of Gravanis et al. are smaller than results given by Fjær et al. At each bottomhole pressure, sand mass given by the model of Gravanis et al. increases until it reaches constant value. On the other hand, the one of Fjær et al. increases continuously because the process of erosion repeats which causes the radius of cavity and production zone to increase, so sand production mass increases continuously. These two models show different results although they use the same physical mechanisms: the fluid flow can be described by Darcy’s law and the driving mechanism for continuous sand production is erosion from plastified material in the vicinity of the production cavity under hydrodynamic forces. However, the model of Fjær et al. uses a step of instantaneous sand production which is not considered in the model of Gravanis et al. Moreover, while Gravanis et al. divide the whole region into a plastic region where Mohr-Coulomb failure criterion is considered, and an elastic region where Hooke’s law is applied, Fjær et al. consider only the plastic region. Fjær et al. simply assume that the stiffness of the rock remains constant until the porosity has reached a critical value ϕcr. At that point, the entire part of the rock that has been producing sand collapses, and the remaining solid material in that part is produced in one burst. As a result, the equations used in these two models as well as the associated experimental parameters are quite different and eventually led to quite different results.Finally, it is important to recall that, due to the lack of validation data in reality, it is impossible to choose which model to use in practical use. At this state, these results can only be used for illustration in study, and to demonstrate the possibility of using Hydromechanical models in real cases to predict eroded sand mass over time. In the future if real data is available, these models can be revised and recalibrated. Unfortunately, until now, measuring sand mass rate is technically impossible in the oilfields. An idea was proposed to solve this problem; the determination of real sand rate may be based on the erosion rate of the choke. However, this study has never been realized and may constitute a new objective in our future study.


## 4.1. Calculate Critical Bottomhole Pressure Using Geomechanical Model of Willson et al. [1]

The critical bottomhole pressure CBHP is important information during production phase because it indicates the lowest bottomhole pressure for sand production to not occur. For a specific well in oil field, the only production data that we control is the bottomhole pressure, which is adjusted using surface chokes. The pressure is controlled; hence the flowrate is controlled. For this reason, before predicting the sand production using Hydromechanical models of Gravanis et al. and Fjær et al., we firstly calculate the CBHP using Geomechanical model of Willson et al.The testing showed that a relationship between the effective in situ strength of the formation U and the TWC (Thick wall cylinder) strength relative to a specimen with an OD/ID ratio (Outer/Inner diameters) would be equivalent to(30)UCS=bf∗1.55∗TWC,where bf the boost factor is often taken as bf=2 for cased–hole perforated wells. Equation (1) considers effect of reservoir pressure decline [9] which can be accounted for by updating in situ stresses, while the vertical stresses are usually kept constant:(31)σh=σh+A·ΔPσH=σH+A·ΔP,where A=α1-2ν/1-ν with α is the Biot’s constant and ν is coefficient of Poisson, ΔP=Pc-Pi with Pc is current reservoir pressure, and Pi is initial reservoir pressure.Figure10 presents the bottomhole pressure Pw versus reservoir pressure Pr. In Figure 10 we also present the CBHP when the UCS changes. For example, for UCS = 1950 psia, if reservoir pressure Pr is 6000 psia, the CBHP is about 1500 psia; if the bottomhole pressure Pw is lower than this value 1500 psia, sand production will occur.Figure 10
Prediction of sand production with variation of UCS (Unconfined Compressive Strength).The results showed that if UCS gets smaller, the sand production window is greater. This can be explained by the reason that the formation rocks around the wellbore are weakened when UCS decreases. Thus, risk of sand production increases.We consider a case study with UCS = 2450 psi,Pr = 4000 psi. From Figure 10, we see that if the current reservoir pressure is 4000 psi, then the critical bottomhole pressure is 2020 psia. For sand production to not occur, the bottomhole pressure must be greater than 2020 psi. Otherwise, to increase the flow rate, we must increase DP and if it is greater than CDP, it means that bottomhole pressure must be lower than critical bottomhole pressure (1000 psi); then the sand will be produced. So there is a compromise between increasing flowrate (production requirement) and sand production.


## 4.2. Calculate Sand Production Mass Using the Models of Gravanis et al. and Fjær et al.

We use Hydromechanical Erosion models of Fjær et al. [5] and Gravanis et al. [6] to calculate sand production mass versus time. However, we do not have experimental data of Well X1 to calibrate empirical parameters, so we must use empirical parameters already calibrated in Section 3. In Figure 11 the results of sand mass produced over time calculated by the two models are presented. We also did a sensitivity study of the effect of drawdown pressure. Drawdown pressure, which is the difference between the reservoir pressure and the bottomhole pressure, is the most important data that we control during production phase. We control the drawdown; consequently, we control the production rate. The drawdown is typically controlled by surface chokes. The more the choke is open, the bigger the bottomhole pressure is and, hence, the smaller the drawdown is, and vice versa.Figure 11
Results of two models in the case study of Well X1.Figure11 shows that sand production mass increases when drawdown pressure increases. However, results given by the model of Gravanis et al. are smaller than results given by Fjær et al. At each bottomhole pressure, sand mass given by the model of Gravanis et al. increases until it reaches constant value. On the other hand, the one of Fjær et al. increases continuously because the process of erosion repeats which causes the radius of cavity and production zone to increase, so sand production mass increases continuously. These two models show different results although they use the same physical mechanisms: the fluid flow can be described by Darcy’s law and the driving mechanism for continuous sand production is erosion from plastified material in the vicinity of the production cavity under hydrodynamic forces. However, the model of Fjær et al. uses a step of instantaneous sand production which is not considered in the model of Gravanis et al. Moreover, while Gravanis et al. divide the whole region into a plastic region where Mohr-Coulomb failure criterion is considered, and an elastic region where Hooke’s law is applied, Fjær et al. consider only the plastic region. Fjær et al. simply assume that the stiffness of the rock remains constant until the porosity has reached a critical value ϕcr. At that point, the entire part of the rock that has been producing sand collapses, and the remaining solid material in that part is produced in one burst. As a result, the equations used in these two models as well as the associated experimental parameters are quite different and eventually led to quite different results.Finally, it is important to recall that, due to the lack of validation data in reality, it is impossible to choose which model to use in practical use. At this state, these results can only be used for illustration in study, and to demonstrate the possibility of using Hydromechanical models in real cases to predict eroded sand mass over time. In the future if real data is available, these models can be revised and recalibrated. Unfortunately, until now, measuring sand mass rate is technically impossible in the oilfields. An idea was proposed to solve this problem; the determination of real sand rate may be based on the erosion rate of the choke. However, this study has never been realized and may constitute a new objective in our future study.


## 5. Conclusion

This study combines Geomechanical model and the two Hydromechanical Erosion models of Fjær et al. [5] and Gravanis et al. [6] aiming to predict sand production. The proposed model can not only predict the critical pressure for onset of sand production but also estimate sand production mass and the variation of porosity in time and space. In addition, the study also gives detailed calculation steps for the two Hydromechanical models. Although each model has assumptions and empirical parameters, we can adjust them if we have experimental data. The results of this study will help to make an effective production planning that can avoid sand production when reservoir pressure declines. Hydromechanical model is still very new and not yet widely used, due to several reasons such as complicated calculation methods and parameters that need to be experimentally adjusted in practice or in the laboratory. Therefore, this initial computational workflow using MATLAB will significantly contribute to the development of future research in this domain in Petroleum Engineering.The experimental data (sand production’s mass in function of time) is rarely done not only in Vietnam but also in the world (only some laboratories and authors mentioned in the paper have done this kind of experience because of their research in this field). Hence, it is currently impossible to obtain empirical data for Cuu Long basin (sand mass in function of time), not only data in laboratory but also production data due to difficulties in collecting and measuring sand production rate. Therefore, application of the calibrated model for Cuu Long basin in this study is destined to demonstrate the possibility of using Hydromechanical models in real cases and the result can be used for illustration in research. In the next studies when we have available data from Cuu Long basin, we can do a revision/recalibration of these models. We already thought about the determination of real sand rate data based on the erosion rate of the choke. This idea will be the objective of our future study.



---

*Source: 2195404-2017-10-17.xml*