Microscopic sand production simulation and visual sanding pattern description in weakly consolidated sandstone reservoirs

Chang-Yin Dong ,Bo Zhou ,Fan-Sheng Huang ,Lei Zhang ,Yi-Zhong Zhao ,Yang Song ,Jun-Yu Deng

a Key Laboratory of Unconventional Oil & Gas Development,Ministry of Education,Qingdao,266580,Shandong,PR China

b School of Petroleum Engineering,China University of Petroleum (East China),Qingdao,266580,Shandong,PR China

c Research Institute of Petroleum Engineering Technology,Shengli Branch,SINOPEC,Dongying,257000,Shandong,PR China

d Research Ltd.of Engineering Technology,CNPC,Tianjin,345000,PR China

Keywords:Weakly consolidated reservoir Particle-scale microstructure model Microcosmic sanding process simulation Visual sanding cavity description Sanding prediction Sand control optimization

ABSTRACT To visually describe the sanding pattern,this study constructs a new particle-scale microstructure model of weakly consolidated formation,and develop the corresponding methodology to simulate the sanding process and predict sand cavity shape.The microstructure model is a particle-objective model,which focuses on the random sedimentation of every sand grain.In the microstructure,every particle has its own size,sphericity and inclination angle.It is used to simulate the actual structure of cemented granular materials,which considers the heterogeneity and randomness of reservoir properties,provides the initial status for subsequent sanding simulation.With the particle detachment criteria,the microscopic simulation of sanding can be visually implemented to investigate the pattern and cavity shapes caused by sand production.The results indicate that sanding always starts initially from the borehole border,and then extends along the weakly consolidated plane,showing obvious characteristic of randomness.Three typical microscopic sanding patterns,concerning pore liquefaction,pseudo wormhole and continuous collapse,are proposed to illustrate the sanding mechanism in weakly consolidated reservoirs.The nonuniformity of sanding performance depends on the heterogeneous distribution of reservoir properties,such as rock strength and particle size.Finally,the three sanding patterns are verified by visually experimental work.The proposed integrated methodology is capable of predicting and describing the sanding cavity shape of an oil well after long-term sanding production,and providing the focus objective of future sand control measure.

1.Introduction

The production of sand (also known sanding) is a universal phenomenon that occurs during the hydrocarbon production in weakly consolidated sandstone reservoirs.Sand production may cause many undesirable problems including,but not limited to,filling and blocking of wellbores,wellbore collapse,plugging of perforations or flowlines,erosion failure of downhole and surface facilities,and additional cost of remedial and cleanup operations(Rahmati et al.,2013).Nowadays,it has been well accepted that certain amount of sand production can be helpful for improving near-wellbore formation conductivity and productivity.But poor management or uncontrol of sand production can lead to problems of sand overproduction with severe consequences aforementioned(Garolera et al.,2019).For decades,the problem of sand production has been approached with the use of sand avoidance techniques,called sand control methods.These methods involve mechanical method to physically block sand grains with screen or gravel packing,and chemical method to reconsolidate the formation and prevent sand from flowing with resin or other chemical agents(Dong et al.,2009;Parlar et al.,2015;Garolera et al.,2019).

As the foundation of sand control design,sand production prediction is extremely important work to clarify the sanding mechanism and law,which can offer strong support for sand control optimization(Dong et al.,2015,2017).Over the decades,due to its importance in hydrocarbon production industry,systematic efforts have been made to develop the numerical methods for sand production prediction.The early studies of sand prediction focused on the macroscopic sand production (Weingarten and Perkins,1995;Nouri et al.,2002;Papamichos and Furui,2019),which mainly depended on the solid failure mechanics of compression,tension and shear effect (Risnes et al.,1982;van Den Hoek et al.,2000).These methods concern the use of rock strength,formation stress distribution and failure criterion of rock to calculate the critical drawdown pressure for onset of sand production(Wan and Wang,2004a;Rutqvist et al.,2012).This type of sand prediction models can only answer under what production conditions (production rate or bottomhole flowing pressure) the wells begin to produce sand,but is unable to predict the sand production rate and volume quantitatively (Crook et al.,2003;Dong et al.,2017).In order to make volumetric predictions of sanding,series of numerical models for sanding simulation and prediction were developed continuously.The early sand production and erosion models for predicting sand production in unconsolidated formation was presented by Vardoulakis et al.(1996).This mathematical model was based on mass balance of the produced solids and flowing fluid in porous media,and can be used to predict the sanding rate with production time.This model was developed furtherly by Papamichos and Stavropoulou (1998) by coupling the poro-mechanical behavior of the solid-fluid system with the erosion behavior of the solid matrix due to fluid flow.Then,further improvements were made by Papamichos et al.(2010),Papamichos and Furui (2019),and Gravanis et al.(2015).This group of models provide mainly an access to predict sand production rate based on flow-coupling and stress analysis.In recent years,the new developed sanding simulation methods,in general,are categorized as continuum and discontinuum approaches.In the continuum approach,cemented granular materials of sandstone formation are treated as continuous matters in deriving the governing differential equations(Morita et al.,1989;Nouri et al.,2006,2007;Papamichos et al.,2010;Hussein et al.,2018).This type of sanding prediction models depends on the presupposition of continuous and homogeneous medium of formation.Garolera et al.(2020) presented a micromechanical approach based on zero-thickness interface elements for modelling advanced localization and cracking states of cemented granular materials.The increase in computational resources due to discretization of the grain area is compensated by modelling the remote area with continuum finite elements.The model is capable of modelling localization of deformation,disintegration and cracking of cemented rock formation.However,the heterogeneity is not fully considered in his model.The typical method of discontinuum approach is discrete element method(DEM) (O'Connor et al.,1997;Wan and Wang,2004b),which is a useful access to simulate sand production,especially to understand the mechanism of sanding(Rahmati et al.,2013;Ghassemi and Pak,2015;Han and Cundall,2017).However,it is difficult to describe the heterogeneity of rock properties and its application is limited in large scale simulation(Climent et al.,2017).To date,similar studies also have been done to sand production prediction for gas hydrate bearing sediments with particular consideration of hydrate decomposition (Cheng et al.,2010;Li et al.,2013;Pinkert and Grozic.,2015;Uchida et al.,2016;Kajiyama et al.,2017;Dong et al.,2019).Generally speaking,the prediction results of abovementioned continuum and discontinuum methods just show the macroscopic sanding law with oil or gas production,but are not capable of describing the sanding pattern and cavity shape of nearwellbore formation and their heterogeneity.

In an actual well in a weakly consolidated sandstone reservoir,after a long-term production with sanding,the continuous removal of solid from formation tends to cause the obvious change of microscopic structure and porosity-permeability properties (Wan and Wang,2004b;Dong et al.,2017,2020).From the engineering point of view,one of the interesting issues the petroleum engineers care is what the sanding pattern and cavity shape would be after a given period of sanding production,the recognition of which can essentially help to optimize sand control design.Unfortunately,although major improvements have been achieved in the past decade,it is still a difficult and challenging issue to describe the sanding pattern and cavity shape of an actual sanding well before sand control design and implementation.What we just know is that the sanding cavity may be very complex and nonhomogeneous depending on the characteristic of formation properties (Dong et al.,2020).

In this work,we construct a new particle-scale microstructure model of a weakly consolidated formation core and the subsequent numerical methodology to simulate the microscopic sanding process during production.The proposed microstructure model describes the heterogeneity of formation using random function,considering the particle size distribution (PSD)curve,well logging data,formation stress distribution and dynamic rock strength.This integrated methodology provides a new access to predict and describe the sanding cavity shape of an oil well after long-term sanding production,which is essentially meaningful for sand control design optimization.

2.Microstructure and characteristic model

2.1.Basic microstructure model

Microstructure model is the foundation for microscopic sanding simulation.Its task is to construct the microscopic structure as precisely as possible,depending limited data and information which can be obtained from oilfield.To approach to the truth of the core,the microstructure should consider PSD of formation sand,rock strength,stress field,fluid properties,etc.The microstructure model in this study is particle-size-scale model.The basic assumptions for the model construction are listed as follows:(1)the model is a 2D model and the particles are represented by plane graphics.(2) Every irregular particle has its own particle size,sphericity and inclination angle,which will be randomly defined in detail later.(3) Every irregular particle has an excircle called subjunctive equivalent circle (SEC).The SECs of actual sand grains are tangent to the adjacent SECs.(4) There is a stress chain and a strength chain between every two particles.

Fig.1a show scanning electron microscope (SEM) image of an actual weakly consolidated core.The sand grains are of different sizes,shapes and promiscuous locations,all of which seem to be unordered.The idealized 2D microstructure model shown in Fig.1b is the model we will construct here.The model is a particleobjective model,which focuses on the random settlement of every sand grain.As seen from Fig.1b,the sand grains with different sizes show obvious randomness,but the general size characteristic statistically accords with the PSD curve.In order to consider the randomness of grain size under the controlling of given size distribution,we use the method of random settlement simulation (Bo et al.,2003) to simulate the final location of every grain.This method employs a container and simulates the process of sand grain settling into the container and final accumulation state of all grains.The first step of the simulation is to pick up a SEC grain random and throw into a container.Then,SEC grain tends to settle at a random location but comply the rule of that it will fall down and reside at the location with stable supporting relationship with other grains present in the containers.Repeating the above process,a basic microstructure can be obtained since the accumulation of particles is random and will eventually form a stable accumulation state.Resultantly,the size and location of every grain can be obtained accurately.What should be noted is that,in Bo's model,the particles are regarded as spherical.In this study,the sand grains are not simplified as circle,but with different inclination angle and sphericity.we assume a subjunctive equivalent circle(SEC) for every actual sand grain,shown by the yellow circle in Fig.1b.SEC is defined as the standard circle surrounding the outer boundary of the grain,and is used to simulate the random piling of the grains.In the actual core shown in Fig.1a,the sand grains contact other grains closely with cement.

Fig.1.Microstructure model of a weakly consolidated sandstone core.

To simulate the random settling of sand grains,one of the preparation works is to define the properties of all sand grains,including size and shape.From oilfield onsite,the information which can reflect the characteristic of sand property is PSD curve.PSD curve provides the grain size distribution with the corresponding weight percent (Fig.2a).Obviously,in the simulation of grain settlement,what we need is the number of the grains with every size range.The following formula can be used to convert PSD data from weight to number with given size range.

wheremis the total mass of sand sample,kg;ρgis the density of sand material,kg/m3;diis the grain size with order ofi,m;wiis the weight fraction of grain size with order ofi,dimensionless;Niis the corresponding numbers of grain with order ofi.

Fig.2 shows a case of converting the PSD curve from weight percent to number percent.The median size of the sand sample is 0.24 mm.The sand sample is from Gudong oilfield located in Dongying,Shandong Province,China.Please note that this PSD curve will be used in the case study later.Fig.2a is the initial PSD curve with weight percent.The converted PSD curve with number percent is shown in Fig.2b.With the same weight percent,the amount of grains with small size is much higher than that of grains with large size.As shown in Fig.2,the location of peak value of bar group moves from about 0.35 mm in Fig.2a to about 0.25 mm in Fig.2b.In the grain settlement simulation,all of the sand grains can be prepared ready according to the amount of total grains and the size distribution shown like Fig.2b.

Except grain size,sphericity and inclination angle (α) are other two key parameters to influence the detachment of grains from the core body.In the actual core shown in Fig.1a,the sphericity and inclination angle seem to distribute randomly.In order to take this feature into consideration and reduce the complexity of the simulation,total 6 types of grain sphericities and 4 types of inclination angles are preset as shown in Fig.3.The preset sphericities involve 1.0,0.9,0.8,0.7,0.6,and 0.5 with the distribution probability (P) of 0.10,0.15,0.25,0.25,0.15,and 0.10,respectively.The series of inclination angles are set as 0°,45°,135°,and 90°with distribution probability of 0.20,0.30,0.30,and 0.20,respectively.

After all of the grain properties are preset ready,the grain settlement simulation can begin similarly to the procedure mentioned above.For one step of simulation,a grain is picked up randomly,with randomly assigned sphericity and inclination angle by preset values,and then“thrown”into the simulation container.As a result,the microstructure model like what is shown in Fig.1b can be obtained.In this model,the location,sphericity and inclination angle are of the feature of randomness,but the statistics of the grain size is accordance to the PSD curve of the sand from target reservoir.

Fig.2.Converting sample of PSD curve from weight percent to particle number percent.

Fig.3.Preset grain sphericities and inclination angles and their distribution probability.

2.2.Intergranular strength distribution

For weakly consolidated sandstone core,the sand grains are cohered together with cement to each other.The intergranular cementation strength is another key property to determine the detachment performance of sand grains.Under normal circumstances,part of the stress in the formation is borne by the pore pressure and part of the stress in the rock mass is borne by the skeleton.The stress borne by the rock skeleton is transferred between the particles in contact with each other through the cementation between the particles.For a certain particle,it is connected to the surrounding particles through a cement and fixed on the solid skeleton.Such a cemented connection is called a connection bond.Each connection bond has a certain connection strength and can transmit pressure and shear force at the same time.In the microstructure model shown in Fig.1b,the next important work is to determine the intergranular strength distribution over all of the grains.

From the field onsite,well logging data is always easy to be acquired for rock strength capturing.According to rock density and acoustic logging data,the cohesion strength,tensile and compression strength of rock can be calculated out by empirical formulas easily(Dong et al.,2017).For a vertical well,the vertical distribution of rock strength along the wellbore axis can be used to determine the vertical intergranular strength in the microstructure model.As shown in Fig.4a,the vertical rock strength distribution from well logging data indicates the heterogeneity of the formation,which is directly tested by well logging.So,the description of intergranular strength heterogeneity along the boundary,i.e.the wellbore,can be achieved successfully.Generally,the maximum horizontal and minimum principal stresses are obtained according to the actual insitu stress test results.Then the Heim model of density logging data is used to calculate the vertical principal stress,and then the horizontal tectonic stress coefficient is obtained through the twophase unequal stress model.The horizontal stress structure coefficient is used to calculate the principal stress profile of the original formation,which is used to predict the in-situ stress of other wells in the same area.

Logging data is the main information for engineers to investigate the properties of reservoir.The fluctuation character of the logging curve,which can be seen easily from logging curves of actual vertical and horizontal wells,indicates the heterogeneity of reservoir rock not only in the vertical direction but also in the horizontal direction.In the case of a vertical well,which is what this study aims at,the logging data can be used to represent the vertical heterogeneity of rock strength.But there is no applicable data from oilfield to indicate the horizontal heterogeneity of reservoir,which is also a key issue for the construction of microstructure model in this study.To solve this problem,we apply the method of random distribution and use vertical logging data to describe the horizontal strength distribution.For the horizontal direction at a certain depth in the vertical well,the intergranular strength and its heterogeneity can be characterized by random function which fits normal distribution.The average value of random strength can be determined as same as the value from the vertical strength distribution curve at the same depth.For the given depthh,the strength of the grain with orderican be calculated out using formula as follows:

whereXsis the non-homogeneous coefficient,dimensionless;Kh0is the average value of strength along the horizontal direction at depthh,MPa;KhminandKhmaxare respectively the minimum and maximum values of strength at depthh,MPa;Kiis the strength of grain with orderi,MPa;Xiis a dimensionless random factor with random value within 0-1;Nis the grain number along the horizontal direction at depthh.

It is assumed that the grain strength along the horizontal direction conforms to a normal distribution.In order to determine the value ofXi,a virtual data box forXito be picked up is shown in Table 1.

Table 1 The virtual data box and Xi distribution.

Fig.4.Logging curve and microstructure model construction.

In the data box shown in Table 1,there isMtypes of value forXi,which isXk1toXkM,respectively.jis the order number of the type of the value.Nkjis the amount of data with value ofXkj.The value with orderjis determined by:

The amount of data with value ofXkjis calculated by:

where σ is the coefficient of normal random distribution,with recommended value of 0.2;f(j) is the normal distribution factor,dimensionless.

In the constructed data box by Table 1 and Eqs.(5)-(7),the total amount of data is justN:

During the construction of grain strength distribution with Eq.(4),a random numberRndshould be generated firstly by computer.The order number of values in Table 1 is calculated by Eq.(9) as bellow:

Keepingjas integer,the valueXkjwith orderjin Table 1 is picked out of the box,meanwhile the amount of the valueXkj,that isNkj,changes toNkj-1.Repeating the above process of picking data until all of the data are used up.

The non-homogeneous coefficient,Xs,is used to express the degree of strength heterogeneity at the horizontal or plane direction of the formation.This coefficient determines the changing range of the random value due to the heterogeneity,which can be calculated by the vertical heterogeneity of the strength,using the following formula:

whereSis the average value of rock strength within one pay zone,MPa;KminandKmaxare the minimum and maximum values of strength of the pay zone,respectively,MPa.

According the PSD curve (Fig.2b) and strength distribution curve (Fig.4a),the resultant microstructure is shown in Fig.4b,which is characterized by the inner cohesive strength.The longitudinal axis represents the vertical depth of the zone and the transverse axis is the distance from the center of the wellbore.The different color indicates various strength distribution.Deep color represents high strength.It can be seen that,along the vertical direction,the strength distribution shows obvious randomness,but its general rule along depth is controlled by the strength distribution shown in Fig.4a.

Fig.5.Random distribution of cohesion strength at different locations.

From the result shown in Fig.4b,the intergranular strength distribution at any depth or radius can be obtained.Fig.5a shows the strength distribution along the vertical direction with a distance (radius) of 0.3 and 0.5 m.Along the vertical direction,the strength distribution shows obvious randomness,but its general rule along depth is controlled by the strength distribution shown in Fig.4a.Fig.5b shows the strength distribution along the horizontal(radial) direction at a depth of 1312.5,1313.8 and 1315.1 m.The strength varies randomly within a certain range defined by Eq.(2)to Eq.(10),but the average value is equal to the value at the strength distribution curve shown in Fig.4a at the same depth.It can be concluded that the microstructure model not only considers sand grain properties and the strength distribution with PSD curve and well logging information,but also considers the randomness and heterogeneity of these properties along the vertical and horizontal directions.It should be stressed again that the model provides important foundation for further model construction and sanding simulation.

2.3.Intergranular stress model

As well known,Intergranular stress is also an important factor influencing grain detachment.In an actual vertical well placed in 3D cylindrical coordinate system,the macroscopic principal stress concerns vertical stress,radial stress and tangential stress,which can be obtained from well logging data using relative stress model(Risnes et al.,1982).For 2D microstructure of this study,the macroscopic principal stress is simplified to vertical stress and radial stress.Here,the tangential stress in 3D model is ignored.

In the 2D microstructure,every sand grain is connected to other grains by cementation and attached to the body of core.The cementation between two grains can be thought as a chain.The chain has a strength and also can transmit stress.The intergranular stress can be decomposed as normal force,Fn,along the straight line crossing the geometric center of two grains,and tangential force,Fs,perpendicular to the above-mentioned straight line(Fig.6).The intergranular stress model is used to acquire the distribution of normal stress and tangential stress of any two grains.For the 2D model,with macroscopic vertical stress,σv,and radial stress,σr,the stress balance equation of boundary grains can be expressed as:

wheremis the amount of grains;is the normal stress of grain with orderi;is the tangential stress of graini.Setting the normal stress and tangential stress of each chain as solving target,the following equations can be obtained:

wherenis the amount of chains.The stress balance equation are as follows:

whereaijis the coefficient to characterize the connection between two grains.When two grains are connected by a chain,aij=1;otherwise,aij=0.is the stress of boundary grain.When the grain is not boundary grain,set the stress of the boundary grainto 0.Finally,Eq.(15) can be expressed as follows:

whereAis non-singular matrix,A-1is the inverse matrix ofA.WhenAis a singular matrix,A-1is the generalized inverse matrix ofA.The intergranular stress distribution can be obtained by solving the equations.

Fig.6.Force decomposition of two grains in the 2D microstructure.

3.Sanding simulation method

3.1.Dynamic strength weakening model

For sandstone reservoir,especially a weakly consolidated sandstone reservoir,many studies have demonstrated that an increase in water saturation tends to reduce the rock strength.The mechanism of rock strength weakening has been explained that water can dissolve and soften the cementation in the sandstone rock,and cause the cementation strength to decrease(Petrobel and Wally,2007).In an actual well in the weakly consolidated sandstone reservoir,as production continues,the water saturation will increase due to the natural migration of oil-water interface.The law and performance of water saturation can be approximately estimated by reservoir numerical simulation.Studies has shown that the rock strength weakening depends on the degree of water saturation improvement and the time of core being soaked in water(Skjaerstein et al.,1997;Dong et al.,2017).

Dong et al.(2017) performed a series of experiments to investigate the performance of sandstone rock strength weakening with water saturation.They put forward an empirical model to describe the relationship of rock strength and water saturation.Based on their studies,an integrated dynamic strength weakening model is developed to predict the rock strength performance for sanding simulation,in which the effect of water on rock strength in sanding simulation is considered.The model is expressed as below:

wheretis the production time,d;Swtis the reservoir water saturation at timet;Sw0is the initial water saturation;tiis the final time that the water saturation increases to,dSw;α and β are fitting coefficients;Kcminis the compressive strength ratio of final state after saturation,fitted from experiment data;Ktis the dimensionless intensity ratio that varies with saturation time;Scis the rock strength,MPa;Sc0is the rock strength under the initial water saturation in the initial stage of production,MPa.

According to the result of Dong et al.(2017),although the strength of sandstone tends to decrease with an increase in water saturation and the soaking time,the strength weakening does not persist for ever.When the strength decreases to a certain threshold value,it will not change again and keep stable.This threshold strength is really the minimum value of strength after water weakening,which is called ultimate strength.The abovementioned coefficientKtminis just the ratio of ultimate strength to the initial strength.This empirical model builds the relationship of rock strength and production time,which makes the sanding simulation method be capable of considering the effect of time on sand production.It should be noted that the relationship between water saturation and production time can be obtained from the overall development plan (ODP) of the reservoir.Normally,in the ODP of a specific reservoir,reservoir numerical simulation can provide the changing performance of water saturation with production time.

3.2.Particle detachment criteria

Particle detachment criteria is used to judge whether a sand grain can detach from the body of core under certain conditions.The detachment of grain from the body depends on so many factors which involve grain size,grain sphericity,angle of inclination,cementation strength,intergranular stress,fluid density and viscosity,in-situ flow velocity and direction,pore pressure and the status of the grain being surrounded by other grains.All of the above properties have been determined in advance by the microstructure model,except the grain surrounding status.Due to the randomness of particle size,sphericity and inclination,the microstructure of sandstone is very complicated.No matter how the particles are attached to the sand body or how the sand particles are surrounded by other particles,it is difficult to answer clearly with a general model.In this study,the simplified structural model summarizes the complicated state where grain is surrounded by other grain to three simple models (case A,case B and case C),as shown in Fig.7.

Fig.7.Structural models of grain detachment.

Fig.7a depicts the first case(case A)of a grain being surrounded by other two grains.This case represents the easiest situation of a grain to detach from the body of the core,where the grain is exposed out of the body with maximum degree of exposure.In this situation of case A,the target grain (with color of blue and red in Fig.7a) is cemented to other grains with minimum intergranular strength due to its very little contact area.Meanwhile,in case A,the target grain will be imposed the maximum fore and torque by the flowing fluid,due to its very high exposure degree value(expressed as symbol ofEin Fig.7).Here,we define exposure degree as the distance from the target grain center to the straight line which connects the center points of the surrounding grains.The schematic drawings of exposure degree under different cases are depicted in Fig.7.In case A,the value ofEis of the maximum value ofEmax,which can be understood as one of the extreme cases.

Another extreme case of grain surrounding is presented in Fig.7c as case C.In this case,the target grain is almost submerged by other grains,which is characterized by the exposure degree of zero.Case C presents the extreme situation in which the target grain may detach from body of the core,but with maximum difficulty.That is because that the target is imposed the maximum cement strength and minimum drag fore by the flowing fluid,meanwhile with the shallowest channel to be produced out.Furtherly in case D shown by Fig.7d,the target grain loses its producing channel with the exposure degree ofE＜ 0,so it is impossible for it to be produced out,even if it can detach from the body of the core.In a word,case D will be judged as no sanding takes place due to the exposure ofE＜0.

Comparing to the two extreme situations of case A and case C,Case B shown by Fig.7b,presents the ordinary situation of grain being surrounded by other grains,which is characterized by the exposure degree ofE＞0 andE＜Emax.Due to its universality and representativity,case B is taken as an example to develop the grain detachment criteria by force analysis.The force of the target grain to hamper its detachment mainly involves the cohesion strength and the possible friction force from the contacted grains.In contrast,the force to induce grain detachment involves the intergranular stress and the drag force from the flowing fluid.According to torque balance,the particle detachment criteria can be expressed as:

whereTiis the torque to induce the grain detachment;This the torque to hamper grain detachment;λ is the coefficient to express the criteria of grain detachment,which is recommended a value of 1.2 here.The detailed calculation models and formulas can be found in early studies(Garolera et al.,2019;Wan et al.,2002)and will not be repeated here.

3.3.Sand production simulation method

With the microstructure model and the defined distribution of all properties,sanding simulation can be performed step by step,the principle and method will be presented in this section.Fig.8 shows the basic principle of sanding process simulation with three steps,where different colors are marked with grains to distinguish different types of grains.

Fig.8a shows the initial state,also the first step of the simulation process.With the assumption of left edge as out boundary,the grains withlight bluecolor are the boundary grains.According to sanding mechanism (Dong et al.,2017;Garolera et al.,2019),the first group of sand grains which fall down from the body come from the boundary grains.Using the grain detachment criteria(Eq.(16)-18),the detaching grains,the grains that will detach from the body,can be identified out,which is colored with red and marked with number 1 in Fig.8a.

After the first group of detaching grains is produced from the body,the boundary grains will change.Some grains ever hidden by the produced grains become new boundary grains,which is colored with color of light blue and red as shown in Fig.8b.Similarly,the new detaching grains will come from the new boundary grains,which can be picked up also using Eq.(16)-18.For easy of distinguishing,the second groups of detaching grains are marked with color of red and number of 2 in Fig.8b.It is unnecessary to go into details that Fig.8c shows the next cycle after the second step.With the above step by cycle,the sanding process simulation can be performed continually until no sanding takes place or the process is stopped manually.Finally,the whole method and model of microstructure generation and sanding simulation have been programmed into a computer module for further case study.The computation module is called scSAND,which is the key module for sand production simulation in the software platform of SandControl Office.The module and platform were self-developed by the team of the authors.

4.Results and discussion

4.1.Microscopic sanding pattern simulation

The oil reservoirs in Gudong oilfield are typical weakly consolidated sandstone reservoirs.The wells in these reservoirs face serious problems of sand production.We select three wells from one of these reservoirs as cases to simulate the sand production process and analyze the sanding cavity performance.Well A is a vertical well,and main basic data used for simulation is shown in Table 2.Please note that the PSD curve of well A is shown in Fig.2a,which is used to explain the basic microstructure model construction in subsection 2.1.The logging and rock properties are presented in Fig.4a,also for microstructure model construction in subsection 2.2.

Fig.8.Diagrammatic sketch of principle of sand production simulation method.

Table 2 The main basic data of well A.

Using the data listed in Table 2,the PSD information shown in Fig.2a and the rock strength curve shown in Fig.4a,the microstructure model of the target zone of well A can be constructed by the method proposed in this work,which is presented in Fig.4b,and the characteristics of which also has been analyzed in subsection 2.2.Based on the microstructure model of well A,the result images of sanding simulation are shown in Fig.9 with production time of 7,65 and 140 d,respectively.

In Fig.9,the area withbluecolor represents the sanding cavity.Other area with color of green to red is the zone with no sanding,the colors in which indicates different rock strength.As can be seen from Fig.9a,at the early stage of production,the detachment of sand grains occurs along the wall of the perforations,and extends around the perforations outward.It is clearly found that the frontier of the cavities formed by detachment of sand grains seems to be very irregular with some feature of randomness.Comparing the sanding cavity of production time of 7,65 and 140 d shown in Fig.9,as production continues,the sanding frontier extends upward,downward or forward,but generally along the radial direction.Looking carefully at the sanding frontier area shown in Fig.9b and c,there are many narrow and irregular extending channels,which look like wormholes in the earth.This sanding pattern implies that from the initial detachment point of grains,the microscopic cavity tends to extend along the profile or path of relatively weak consolidation,where the disaggregation of cemented grains takes place easily.Apparently,the primal cause of the irregular extension performance is the heterogeneity and randomness of rock properties defined be the previous microstructure model.In this study,we call this type of sanding pattern as pseudo wormhole sanding(PWS)pattern,the detailed description of which will be illustrated in next subsection 4.4.

In the area close to the wellbore and perforations,detachment of sand grains seems to be very serious.With the small wormholes extending and the connecting of different wormholes,relatively big cavity comes into being.This sanding pattern looks like formation collapse very much.Here,we name this type of sanding pattern as continuous collapse sanding (CCS) pattern.The CCS pattern is characterized by the instant detachment and collapse of grains from rock body on a relatively large profile.Also,the detailed description of CCS pattern will be described in subsection 4.4.As an overall trend,the sanding pattern of well A can be thought as the combination of PWS pattern and CCS pattern.

Fig.9.Visual simulation image of sanding cavity extension of well A.

Another important concern of sanding simulation result is to investigate the longitudinal characteristic of sand cavity along the wellbore axis.In Fig.9,it can be seen obviously that the sanding frontier at deferent depths is strongly jagged,but not parallel.In other words,the sanding frontier show strong nonuniformity along the vertical axis.The general tendency is that the sanding frontier extension of upper section is slower than that of the lower interval.The section with a depth of 1312.0-1312.5 m shows the weakest sand production performance,and the lower section with a depth of 1314.4-1314.8 m shows the worst sanding performance.What is the reason for this tendency?We can find answer from the logging and strength distribution shown in Fig.4a that the strength of the upper section is obvious less than the lower interval.This accordance of initial logging data and the simulation result of sand cavity implies the validity of this model and method,which successfully takes the randomness and heterogeneity of reservoir properties in account at the same time.

4.2.Sanding cavity shape analysis

Except for investigation of microscopic sanding pattern,the more important application value of this presented sanding simulation methodology lies in its capacity to simulate and predict the macroscopic sanding pattern and cavity shape of an actual well,after the technical processing of scale upgrades.

Well B is also a cased-perforated vertical well in Gudong oilfield from the same reservoir as well A,which we chose for further case study of sanding cavity prediction.The target zone of well B covers the depth from 1317.5 to 1320.0 m,with a perforation length of 0.45 m and a shot density of 16 shots per meter.The used PSD information is the same as well A shown in Fig.2a.According to the necessary basic data,the sanding simulation is performed on well B,and the results of sand cavity extension are presented in Fig.9 with different production time of 5,20,40 and 125 d.

From the simulation results shown in Fig.10,it can be observed that well B also displays strongly heterogeneous sanding profile along wellbore.The sanding cavity extension shows characteristic of irregularity and randomness.Relatively bigger wormhole propagation can also be seen as PWS pattern.In fact,it also can be understood as CCS pattern,because the judgement of sanding pattern partly depends on the scale of view one uses.Generally speaking,the upper section close to the top boundary shows the strongest sanding tendency.The section with a depth of about 1319.0 m is another area producing sand seriously.

The macroscopic description of sanding cavity shape is very useful for the following sand control measure decision.For example,one of the principle or rules is that,for an actual sanding well,the section with severe sanding tendency or big sanding cavity should be packed with gravel preferentially to ensure good sand control effect.In order to present the application of the proposed methodology in macroscopic sanding cavity prediction and sand control optimization,another well,well C,from Gudong oilfield is chosen as new case analysis.The thickness of production zone of well C is 10 m,much longer than that of well A and well B for showing macroscopic sand cavity shape clearly.Using the proposed method and self-developed program,the sanding cavity shape of well C is shown in Fig.10a,seeing the area with yellow color.Based on existence of sanding cavity,the gravel-packing process was also simulated visually,the results of which are shown in Fig.10a,b and c with different packing stages.Here,the packed area is marked with the color of red-yellow.What needs to be noted is that the gravel-packing simulation method is not the topic of this study.The objective of its presence here is just to help to illustrate the usage of sanding cavity prediction in sand control optimization.

Fig.10.Visual simulation image of sanding cavity of well B.

Fig.11.Visual simulation image of sanding cavity and gravel packing of well C.

Fig.11a depicts the general sanding cavity of well C.The sanding pattern also is the combination of PWS pattern and CCS pattern.It appears that almost all of the production zone tends to produce sand during oil and water production.Particularly,close to the bottom boundary of the zone,there is a big and deep sanding cavity.As previously mentioned,this position of big cavity will be the focus objective of the following sand control measure.Figs.11b and c depict the medium and final stages of gravel-packing in sand control job without special optimization.It can be seen that at the position of the deep cavity,there will leave an unpacked cavity after gravel packing job (Fig.11c).When the well is re-put into production,the unpacked zone tends to cause instability of packing gravel and further unwanted problem.To fully pack the sanding cavity,the particular optimization under the instruction of sanding cavity prediction may be needed,which is not discussed here.

It should be noted that the simulation result is based on 2D model.Some calculation has been converted from 3D situation of actual well.But,the converting calculation has been simplified greatly,which leads to some deviation obviously.Even so,the results of this work still have essential guidance for sand production prediction and sand control.

4.3.Sanding frontier performance and sensitivity analysis

The sanding simulation method can provide quantitative sand prediction results.The simulation results of the sand production frontier and equivalent sanding cylinder radius varying with production time are shown in Fig.12a and Fig.12b for well A and well B,respectively.Here,the sand production frontier is defined as the maximum distance to the wellbore center of sand grains detachment at a certain time.It can be seen from Figs.12a and b that,in general,the sanding frontier extends outward from the wellbore with production time.It reveals that the formation zone producing sand increases with production continuing with no sand control measure.It can also be seen that,sometimes,the sanding frontier keep no change within a certain period of time,such as from 80 to 120 d in Fig.12a.This interesting performance can be illustrated that the detachment of sand grains does not always propagate forward along the radial direction.Due to the randomness and heterogeneity of the rock properties,the sanding frontier may extend upward,downward,even backward sometimes.

It is clear that the index of sanding frontier can only express the front edge of sanding zone in the formation,but cannot exhibit the severity or the rate of sand production.This is because the value of sanding frontier is really no relative to the sanding pattern and cavity volume within the scope confined by the frontier.Furthermore,sand grain detachment and production tend to propagate along the irregular weak cementation plane,but not a stable and translational profile.

Fig.12.Obtained performance of sanding frontier and equivalent radius with production time.

To express the sanding volume or rate,another index,equivalent sanding cylinder(ESC)radius,is proposed.The ESC radius refers to a special commensurate radius of a cylinder around the wellbore.The volume of the total produced sand is equal to the volume of this cylinder with elimination of the wellbore volume.The performance of ESC radius with production time is also shown in Fig.12.As can be seen in Fig.12a and b,the ESC radius at a certain time is far less than the sanding frontier at the same time.The performance of ESC radius can show sanding severity.For example,the sand production of well A shown in Fig.12a is a little severer than that of well B shown in Fig.12b.It should be note that the above results are from the 2D simulation.In an actual well,the sanding takes place and the frontier extends within a 3D space,so the sanding performance will be much more complicated than that of 2D simulation.

To investigate the effects of geological properties and production conditions on sanding performance,we conducted sensitivity analysis of sand production in well A,using different average cohesion strengths and average liquid production rates.In the sensitivity analysis of average cohesion strength,the values were set as 0.25,0.33,0.46,0.63,and 0.82 MPa.Correspondingly,to keep the consistency of cohesion strength to its vertical distribution curve,the curve shown in Fig.4a is moved parallelly left or right to meet the average value.Setting production time of 100 d,the simulation result curves of the final sanding frontier and ESC radius varying with average cohesion strength are presented in Fig.13a.Similarly,with the basic data of well A shown in Table 2,the production rate was set as 40,45,50,55 and 60 m3/d,respectively,to simulate and analyze the effect of production rate on sanding performance.The corresponding sensitivity analysis curves are shown in Fig.13b.

It can be seen from Fig.13a that sanding frontier and ESC radius decrease with increasing of average cohesion strength.When average cohesion strength increased from 0.25 to 0.82 MPa,sanding frontier experienced a significant drop from 2.39 to 1.09 m.Correspondingly,ESC radius is also slightly reduced.Higher cohesion strength of sand grains in the rock indicates the much higher difficulty for the grain to desquamate from the rock body,which is reflected by sanding rate and sanding frontier or radius.It also can be seen from Fig.13b,the sanding frontier increases with the production rate increasing.When average production rate increased from 40 to 60 m3/d,sanding frontier rose from 1.95 to 2.54 m and ESC radius also increased.High production rate tends to lead to serve sand production due to the high drag force imposed by the fluid on the sand grains.Also,for given reservoir properties,the change of production rate indicates corresponding change of bottomhole flowing pressure or pressure drawdown for production,the relationship of which is determined by the inflow performance relationship(IPR).High production rate will lead to low bottomhole flowing pressure and cause further changes of formation stress distribution,which affect the sanding performance.

4.4.Microscopic sanding patterns and mechanisms

Fig.13.Sanding frontier and ESC radius changes with influencing factors.

Fig.14.Schematic drawing of pore liquification sanding (PLS) pattern.

For deepen the understanding of sand mechanism in weakly consolidated reservoir,according to the phenomenon and sanding pattern obtained from the results of numerical simulations,we summarize three typical sanding patterns concerning pore liquification sanding (PLS) pattern,PWS pattern and CCS pattern.The process and mechanism of these three patterns are reexplained using Figs.14-16,respectively as bellow.It is noted again that PWS and CCS patterns have been observed in the numerical simulation in subsection 4.2.

Fig.14 depicts the sanding pattern and mechanism of PLS pattern.This pattern of sand production occurs in the reservoir with a certain strength,which we call medium-strong cemented reservoir.In this type of reservoir,due to its relatively high cement strength,sand production will be slight and the removal of sand grains from formation will not change the skeleton pattern substantially,but just improve the porosity of the formation near wellbore.Fig.14a shows the initial state of the core.During production,some small grains are easy to detach from the body of the core and are carried by fluid flowing out of the formation.Due to the removal of small grains,the internal space ever occupied by the produced grains become new void space (pore) (Fig.14b),which can be expressed as pore liquification sanding (PLS) pattern.In an actual well,sand production with PLS pattern just changes the porosity and permeability within the formation where sand production occurs (Fig.14c),but does not cause obvious formation deformation.Note that,PLS pattern is difficult to be visually observed,but it can really be demonstrated by experimental simulation later in this study.

The sanding pattern and mechanism of PWS is presented by Fig.15.PWS pattern is prone to take place in the reservoir with weak strength,which we can call weakly consolidated reservoir(here the meaning different from the same term in the topic of this paper).In this type of formation,it is relatively easy to achieve the sand production criteria.Except the pattern of PLS,the more obvious characteristics of the sand production is the PWS pattern.Firstly,along the boundary of core body,which generally refers to the wall of the open wellbore or the perforation channel of an actual well,sand grains detach from the body at the most weakly cemented point and form initial sand cavity.Then,the detachment of grains continues to occur along the weakly cemented plane,meanwhile the sanding cavity extends along the same path.Due to the randomness and heterogeneity of the cementation,as aforementioned in the microstructure model construction,the sanding cavity shows irregular shape.It is very like the hole of worm in the earth,called wormhole as well known.The initial state of weakly cemented core is shown in Fig.15a and the wormhole formed by sand production is shown in Fig.15b.To distinguish the term of“wormhole”used in cold production of a heavy oil reservoir,we call this pattern pseudo-wormhole sanding(PWS)pattern.When PWS pattern is applicable to a sanding reservoir,it indicates that the sand production will change indeed the structure of the formation,as shown in Fig.15c.In this situation,there will be some irregular hole with relatively high porosity and permeability within the formation where sanding.These holes may act as the main channel for the formation fluid to flow and production.Similar studies and analysis about sanding wormhole propagation was also done by Wan and Wang (2004b) using finite element method,which demonstrates the rationality of this sanding pattern.But Wan's results show relative regular extension of wormhole due to its less considering of property randomness and heterogeneity.

Fig.15.Schematic drawing of pseudo wormhole sanding (PWS) pattern.

Fig.16.Schematic drawing of continuous collapse sanding (CCS) pattern.

Another sanding pattern,the CCS pattern,is presented in Fig.16.This type of sand production occurs in the reservoir with extremely low strength,such as quicksand or semi-quicksand formation.There is almost no cementation between the sand grains,also expressed as unconsolidated.In this situation,it is very easy for the grains to reach flowing state.When this type of reservoir produces fluid,the sand will be produced at the same time.From the aspect of mechanism,the sand grains collapse one after one constantly and large cavity tends to form finally,as shown in Fig.16b and c.CCS pattern is really the extreme case of sand production and seem to be seldom in the field onsite,but it is really existing.For example,a few of pay zones of Gudong oilfield,and gas hydrate reservoirs in Shenhu(Cheng et al.,2010;Dong et al.,2019),are thought as semiquicksand or quicksand formation.Certainly,CCS sand production is a disaster for normal hydrocarbon production.Similarly,Garolera's(2019)simulation study shows the sand cavity extends in a very small range around a perforation hole,which is similar to the pattern of CCS.

What calls for special attention is that the patterns of PLS,PWS and CCS proposed in this study are just relative concepts based on the scale of size.This can be illustrated in other words that what pattern a sanding performance should belong to partly depends on the scale we investigate.For example,a sand wormhole with a diameter of 10 mm is certainly categorized as PWS patterns on the scale of eye inspection.But,when observing with much smaller scale,such as enlarging by 100 times,it also looks like wall collapse and can be categorized as CCS patterns.On the other hand,the continuous development of pore liquification will communicate the new created pores together and shape wormhole,which indicates the transition from PLS to PWS pattern.Similar transition of sanding pattern is possible to take place from PWS to CCS.Instant connecting of wormhole will form continuous cavity,impressed by CCS pattern.This transition from PWS to CCS can also be observed in the sanding image of case A,B and C shown in Figs.9-11.

5.Experimental verification

In order to verify the sanding pattern proposed from numerical simulation results,we performed series of experimental simulations of sand production using artificial cement cores with very low strength.The main part of the experimental apparatus is a transparent sanding simulation unit(SSU),the shape of which concerns two types of rectangle and trapezoidal.The size of rectangle SSU is 100 mm (length) × 20 mm (width) × 5 mm,and the size of trapezoidal SSU is 100 mm(length)×15/80 mm(width of narrow/wide section) × 5 mm (thickness).The cover plate of SSU is made of transparent material for capturing image conveniently.Other components of the apparatus involve a constant-flux pump,a fluid tank,a piston container,a camera system,a passed-sand collector,sensors for sensing flowing rate and pressure,flowline and valves.

The principle of the visual sanding experiment is that the constant-flux pump pumps fluid to flow through the prepacked sand media in the SSU to simulate the process of sand production.Before displacement,the sand is mixed with consolidation agent and then packed into the SSU under 80°C for 6 h,a weakly consolidated core is formed after being placed to cool.The strength of the core is controlled by adjusting the concentration and volume of the consolidation agent.When the core is ready,the SSU is connected to the flowline to begin displacement.The fluid is pumped into the SSU and flows through the core and tends to drive the sand grains to move.If the sand detachment criteria are met,some grains with relatively low cement strength will detach from the body of the core and flow out of the SSU with fluid.Continuous sand removal from the core will lead to some types of vacancy,such as continuous cavity,wormhole channel or new-created pore as represented previously.The sanding pattern and cavity shape can be visually observed with the convenience of transparent cover on SSU.For clear investigation,the camera system is used to capture the image of sanding cavity and the video of its extension.

The objective of experiments is to verify the existence of the proposed three sanding patterns,which are summarized in subsection 4.4.The used materials mainly involve fluid,sand and chemical consolidation agent.The displacing fluid is clean water for the convenience of clear observation.The packed sand is manually made by mixing commercial quartz sand with different sizes,according to the PSD data shown in Fig.2a.The chemical consolidation agent is epoxy resin with concentration of 0.1%,0.3%,0.5%,and 0.7% to achieve different cement strengths.

The first series of experiments were performed using rectangle SSU with different concentrations of epoxy resin.Fig.17a presents the resultant picture of sanding pattern with consolidation agent concentration of 0.1%.It can be seen that almost all sand grains detached from the body of the core and were produced out of the unit.Obviously,this sanding pattern can be classified as CCS pattern.When the agent concentration was improved to 0.3% and 0.5%,indicating the core strength also were improved,the final sanding patterns were shown in Fig.17b and c,respectively.It appears that the sand detachment took place and extended along a certain track,which should be the weak cementation plane,to form a channel like an irregular wormhole.It is clear that this sanding pattern can be explained as PWS pattern.When using consolidation agent concentration of 0.7%,as shown in Fig.17d,there was almost no change can be observed on the surface of the core,comparing to the phenomenon with concentration of less than 0.5%.But,in the passed-sand collector,some separate sand grains could be observed clearly.It can be demonstrated that some grains really be produced from the core even if nothing can be seen on the surface of the core.According to the proposed patterns,this sanding pattern belongs to PLS pattern.It might also be noted that the amount of this observed sand is considerably less than that of the experiments using concentration of 0.5% and 0.3%.

Fig.18 shows the resultant pictures of experimental sanding simulations with trapezoidal SSU.The obvious sanding pattern of PWS(Fig.18a) and CCS(Fig.18b and c)were also observed clearly.The zone surrounding with red line is the cavity caused by continuous collapse of grains.During the experiments with trapezoidal SSU,the fluid was injected into the unit from the wide side and flowed out from the narrow side.The arrangement is to simulate the converging flow of near-wellbore formation in an actual well,where the flowing velocity increases with the fluid flows closely to the wellbore due to the shrinking flowing area.

Fig.19 shows the enlarged microscopic image of wormhole forming and propagation captured by the camera system.After the wormhole formed,continuous displacement of fluid will make the wormhole extend along the weakly consolidated plane,which may be somewhat irregular.So,the sanding wormhole normally shown random and irregular shape.In general,from the above experiments,three proposed typical sanding patterns of PLS,PWS and CCS are all be successfully observed with different experiment conditions.This demonstrates that the three sanding patterns do truly exist,and the results of numerical sanding simulation presented by this study are reasonable.

Fig.17.Photographs of different sanding patterns from experimental results with rectangle SSU.

Fig.18.Photographs of different sanding patterns from experimental results with trapezoidal SSU.

6.Conclusions

Sand production in a weakly consolidated sandstone reservoir is a complex phenomenon.The sanding performance depends on various factors involving sand grain size,grain sphericity and inclination angle,intergranular cementation strength,stress distribution,fluid properties,fluid flow field,etc.We present a microstructure model in this work as an effective access to reconstruct the structure of formation near wellbore,which can reproduce the complex behavior of intergranular aggregation and offer the foundation for subsequent sanding simulation.The proposed microstructure model not only fully utilizes the limited information such as PSD curve,well logging data and principal stress,but also specifically takes the heterogeneity and randomness of reservoir properties in account.

Fig.19.Microcosmic image of PWS pattern from experimental result.

Based on the microstructure model,the proposed sanding simulation approach is capable of simulating the microscopic sand production process,depicting sanding cavity characteristic and predicting sanding law quantitatively.In particular,the simulation method and programs can be implemented to visually present the whole process of sand production and sanding cavity propagation.From the simulation results of series of case analyses,the sanding performance shows some interesting and special characteristics corresponding to the heterogeneity and randomness of reservoir properties characterized by the microstructure model.Sensitivity analysis results show that sanding frontier and ESC radius decrease with increasing average cohesion strength and increase with increasing average liquid production rate.For weakly consolidated sandstone reservoir with different properties,it is found that three typical sanding patterns involves PLS,PWS and CCS patterns.In an actual well,the sanding pattern may be the combination of two or three of the typical patterns.The present simulation method and program can be used to predict the sanding cavity shape along the well axis and identify the production profile with severe sanding.This will help to determine the focus target for sand control implementation and improve its effectiveness considerably.

Acknowledgements

This work was financially supported by the National Natural Science Foundation of China (Grant No.51774307,52074331,42002182) and partially supported by Major Special Projects of CNPC,China (ZD2019-184).The first author would like to thank Gudong Research Institute of Oil Production Engineering,Shengli Oilfield,SINOPEC,for providing basic data.